解微分方程的时候,大多数方程都是隐式方法比较稳定。虽然隐式方法需要迭代,写起来比较麻烦,我还是建议尽量使用隐式方法。Adams隐式是一种精度高,稳定性高的算法,属于隐式四阶龙格库塔法的一个特例。我之前写过用python解微分方程的code,这里改成fortran版本
对微分方程dy/dx=f(x),
Adams法,公式
求解的时候用Adams公式构造隐式方程,将y_{n+1}移到右边,然后用牛顿迭代对每个点求解y_{n+1}。下面直接上代码
program potential
implicit real*8 (a-h, o-z)
h = 0.01
icells = 10
y0 = 3.0
y1 = 3.0
y2 = 3.
x0 = 0.0
x1 = x0+h
x2 = x1+h
yi = 0.0 !设置初始值
call RK(y0, x0, h, y1)
call RK(y1, x1, h, y2) !龙格库塔法计算前两个点
do i = 3, 100
xi = i*h
call newton(h, y2, y1, y0, x2, x1, x0, xi, yi) !牛顿截弦法计迭代算后一个点
y0 = y1
y1 = y2
y2 = yi
x0 = x1
x1 = x2
x2 = xi !重新赋值
yyy = 3./(1.+xi**3)
write(*, *)xi, yi, yyy
enddo
stop
end
subroutine newton(h, y2, y1, y0, x2, x1, x0, xi, yi)
implicit real*8 (a-h, o-z)
real*8, external :: fadams
eps = 1.d-15 !精度
i = 0
dx0 = y1/1.d3 !用来计算差商
do while (dabs(dx0)>eps .and. (i<20))
dfx0 = (fadams( yi+dx0, xi, y2, x2, y1, x1, y0, x0, h)-
* adams(yi, xi, y2, x2, y1, x1, y0, x0, h))/dx0
yi1 = yi-fadams(yi, xi, y2, x2, y1, x1, y0, x0, h)/dfx0
dx0 = yi1-yi
yi = yi1
i = i+1
end do
return
end
function fadams(yi, xi, y2, x2, y1, x1, y0, x0, h) !Adams 公式
implicit real*8 (a-h, o-z)
real*8, external :: fun
fn3 = fun(xi, yi)
fn2 = fun(x2, y2)
fn1 = fun(x1, y1)
fn0 = fun(x0, y0)
fadams = yi-y2-h/24.*(9.*fn3+19.*fn2-5.*fn1+fn0)
return
end
subroutine RK(y0, x0, h, y1) !龙格库塔法
implicit real*8 (a-h, o-z)
real*8, external :: fun
rk1 = fun(x0, y0)
rk2 = fun(x0+h/2., y0+h/2.*rk1)
rk3 = fun(x0+h/2., y0+h/2.*rk2)
rk4 = fun(x0+h, y0+h*rk3)
y1 = y0+h/6.*(rk1+2.*rk2+2.*rk3+rk4)
return
end
function fun(x0, y0) !要计算的微分方程右边
implicit real*8 (a-h, o-z)
fun = -x0*x0*y0*y0
return
end