!分别用Euler,Taylor,后向Euler,梯形方法的常微分方程解
program mainimplicit noneinteger(8)::n,m,i,kreal(8)::h,x(100),y1(100),y2(100),y3(100),y4(100),y(20,100),y5(100)real(8),external::ff,f2open(unit=100,file="xable.csv")open(unit=101,file="Euler.csv")open(unit=102,file="Taylor.csv")open(unit=103,file="Back Euler.csv")open(unit=104,file="Trapeziod.csv")open(unit=105,file="Exact.csv")h=0.2n=20m=10                                                                                               !!文件名和变量的定义及赋值do i=1,n+1,1x(i)=(i-1)*hwrite(100,*)x(i)end do                                                                                             !!取n+1个数值点,写入csvclose(100)!!This is the Euler algorithmy1(1)=1write(101,*)y1(1)do i=1,n,1y1(i+1)=y1(i)+h*ff(x(i),y1(i))write(101,*)y1(i+1)                                                                             !!将Euler得到数据导入csvend doclose(101)!! This is the Taylor series methody2(1)=1write(102,*)y2(1)do i=1,n,1y2(i+1)=y2(i)+h*ff(x(i),y2(i))+1/2*h*h*(ff(x(i),y2(i))*(-x(i))-y2(i))write(102,*)y2(i+1)                                                                              !!将Taylor得到数据导入csvend doclose(102)!!This is the back Euler algorithmy(1,1)=1do i=1,n,1y(1,i+1)=y(1,i)+h*ff(x(i),y(1,i))end dodo i=1,n,1do k=1,m,1y(k+1,1)=1y(k+1,i+1)=y(1,i)+h*ff(x(i+1),y(k,i+1))end doend dodo i=1,n+1,1y3(i)=y(m+1,i)write(103,*)y3(i)                                                                            !!将Back Euler得到数据导入csvend doclose(103)!! This is the Trapeziod algorithmy(1,1)=1do i=1,n,1y(1,i+1)=y(1,i)+h*ff(x(i),y(1,i))+h*h*(ff(x(i),y(1,i))*(-x(i))-y(1,i))/2end dodo i=1,n,1do k=1,m,1y(k+1,1)=1y(k+1,i+1)=y(1,i)+h/2*(ff(x(i),y(1,i))+ff(x(i+1),y(k,i+1)))end doend dodo i=1,n+1,1y4(i)=y(m+1,i)write(104,*)y4(i)                                                                            !!将Trapeziod得到数据导入csvend doclose(104)!! This is the exact solutiondo i=1,n+1,1y5(i)=f2(x(i))write(105,*)y5(i)end do
end program main!!dy/dx的function函数
function ff(x,y)implicit nonereal(8)::x,yreal(8)::ffff=-y*xreturn
end
!!exact solution 的function函数
function f2(x)implicit nonereal(8)::xreal(8)::f2f2=dexp(-x*x/2)return
end

学习origin绘制基本点线图对比结果：