2.3多步法求解常微分方程
program main
implicit none
integer(8)::n,m,i,k
real(8),external::fa,fb
real(8)::h,x(200),ya(200),yb(200),yc(200,200),yd(200)
n=100
m=100
h=4.0/n
open(unit=99,file="xable.csv")
open(unit=100,file="Adams-two.csv")
open(unit=101,file="Adams-four.csv")
open(unit=102,file="Adams-Moulton two.csv")
open(unit=103,file="exact.csv")
do i=1,n+1,1
x(i)=(i-1)*h
write(99,*)x(i)
end do
close(99)
!The Adams two-step algorithm
ya(1)=1
ya(2)=1-h*h/2.0
write(100,*)ya(1)
write(100,*)ya(2)
do i=2,n,1
ya(i+1)=ya(i)+h*(3.0/2.0*fa(x(i-1),ya(i))-fa(x(i-1),ya(i-1))/2.0)
write(100,*)ya(i+1)
end do
close(100)
!The Adams four-step algorithm
yb(1)=1
write(101,*)yb(1)
do i=1,3,1
yb(i+1)=yb(i)+h*fa(x(i),yb(i))
write(101,*)yb(i+1)
end do
do i=4,n,1
yb(i+1)=yb(i)+h/24.0*(55.0*fa(x(i),yb(i))-59.0*fa(x(i-1),yb(i-1))&
+37.0*fa(x(i-2),yb(i-2))-9.0*fa(x(i-3),yb(i-3)))
write(101,*)yb(i+1)
end do
close(101)
!The Adams-Moulton two-step algorithm
yc(1,1)=1
yc(1,2)=1-h*h/2.0
do i=2,n,1
yc(1,i+1)=yc(1,i)+h*(3.0/2.0*fa(x(i),yc(1,i))&
-fa(x(i-1),yc(1,i-1))/2.0)
end do
do i=2,n,1
do k=1,m,1
yc(k+1,1)=1
yc(k+1,2)=yc(1,2)
yc(k+1,i+1)=yc(1,i)+h/12.0*(5.0*fa(x(i+1),yc(k,i+1))&
+8.0*fa(x(i),yc(k,i))-9.0*fa(x(i-1),yc(k,i-1)))
end do
end do
do i=1,n+1,1
write(102,*)yc(m+1,i)
end do
close(102)
!The exact solution
do i=1,n+1,1
yd(i)=fb(x(i))
write(103,*)yd(i)
end do
close(103)
end program main
!!exact solution 的function函数
function fb(x)
implicit none
real(8)::x
real(8)::fb
fb=dexp(-x*x/2.0)
return
end
!!-x*y的function函数
function fa(x,y)
implicit none
real(8)::x,y
real(8)::fa
fa=-x*y
return
end