常微分方程数值解法(附录).docVIP

第5章 附录 5.2.1 Euler法 例题5.2.2计算程序 /* demo_euler.c for euler method */ #include math.h double euler(x,y,h,f) double x,y,h; double(*f)(); {y+=h*(*f)(x,y); return(y); } double f(double x,double y) { return(y-2*x/y); } main() { int i; double x, y, h; x = 0; y = 1; h=0.1; for(i=0; i10; i++) {y=euler(x,y,h,f); x+=h; printf(\t%6.2f %12.6f\n,x,y); } return(0); } ============================================== 5.2.2 Runge-Kutta方法 例题5.2.3计算程序 ! runge-kutta method of order 4 ! y = x*y^(1/2) ! y0 = 4 ! external f program demo_runge_kutta implicit none integer,parameter :: n = 100 integer :: i real(8):: y(0:n),x(0:n),h open(2,file=EXA_5_3.txt) x(0) = 2.0; x(n) = 3.0 y(0)= 4.0 h = (x(n)-x(0))/n call rk4(x,y,h,n) write(2,(2f14.8)) (x(i),y(i),i=0,n) end function f(x,y) real(8) :: x,y,f f=x*sqrt(y) end subroutine rk4(x,y,h,n) integer :: i real(8) :: f,h,h2,k1,k2,k3,k4,x(0:n),y(0:n) h2=0.5*h do i=0,n-1 k1=h*f(x(i),y(i)) k2=h*f(x(i)+h2,y(i)+0.5*k1) k3=h*f(x(i)+h2,y(i)+0.5*k2) k4=h*f(x(i)+h,y(i)+k3) y(i+1)=y(i)+(k1+k2+k2+k3+k3+k4)/6.0 x(i+1)=x(0)+h*real(i) end do end ================================================ 5.2.3 微分方程组与高阶微分方程 例题5.2.4计算程序 !!!lorentz.f90!!! program lorentz implicit none real(8) wc,h,z(6,1000),y(6), d(6),t integer i,j,n,m open(2,file=lorentz.txt) y = (/1.,1.,0.,1.,0.,1./) ! 带电粒子初始位置和速度 n = 500 ! 时间步数 m = 6 ! 微分方程维数 h = 0.1 ! 时间步长 t = 0.0 call rk4(t,y,m,h,n,z,d) write(2,(3f14.8)) (z(1,i),z(2,i),z(3,i),i=1,n) end program lorentz subroutine f(t,y,m,d) real(8) y(m),d(m),wc,t wc=1.0 ! 回旋频率取 d(1) = y(4) d(2) = y(5) d(3) = y(6) d(4) = wc*y(5) d(5) =-wc*y(4) d(6) = 0. end subroutine rk4(t,y,m,h,n,z,d) real(8) y(m),d(m),z(m,n),a(4),b(m),t,h,xx,tt a(1)=h/2.0 a(2)=a(1) a(3)=h a(4)=h do 5 i=1,m 5 z(i,1)=y(i) xx=t do 100 j=2,n call f(t,y,m,d) do 10 i=1,m 10 b(i)=y(i)