runge-kutta.txt

两个求解微分方程组 的龙格库塔法程序两个程序可用

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!              Runge-Kutta算法           !!!!!!!
!!!!                     　　　　　　　　　!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      　　SUBROUTINE RKT(y,dydx,n,x,h,yout,F)
DIMENSION Y(n),dydx(n),yout(n),yt(n),dyt(n),dym(n)
real hh,h6,xh,x,h
integer i,n

hh=h/2
h6=h/6
xh=x+hh

do i=1,n
yt(i)=y(i)+hh*dydx(i)
end do
call F(xh,yt,dyt)

do i=1,n
      yt(i)=y(i)+hh*dyt(i)
end do
call F(xh,yt,dym)
do i=1,n
yt(i)=y(i)+h*dym(i)
dym(i)=dyt(i)+dym(i)
end do
call F(x+h,yt,dyt)
do i=1,n
yout=y(i)+h6*(dydx(i)+dyt(i)+2*dym(i))
end do 
end SUBROUTINE RKT
!!!!!!!!!!!!!!!!!!!!!!!
!!!!      main   !!!!            
!!!!               !!!!
!!!!!!!!!!!!!!!!!!!!!!!
program d14r1
      external F
      parameter(n=2)
dimension y(n),dydx(n),yout(n)
x=0.0
y(1)=1
y(2)=1

dydx(1)=y(1)
dydx(2)=-y(2)

do i=1,11
h=0.01*(i-1)
call RKT(y,dydx,n,x,h,yout,F)
write(*,*)'t=',h
write(*,*)(yout(j),j=1,2)
end do 
end program
subroutine F(x,y,dydx)
dimension y(2),dydx(2)
dydx(1)=y(1)
dydx(2)=-y(2)
end subroutine F












　/**
　　***四阶Runge-Kutta法***
　　 
　　 经典格式：
　　y(n+1) = y(n) + h/6 ( K1 + 2*K2 + 2*K3 + K4 )
　　K1 = f( x(n) , y(n) )
　　K2 = f( x(n+1/2) , y(n) + h/2*K1 )
　　K3 = f( x(n+1/2) , y(n) + h/2*K2 )
　　K4 = f( x(n+1) , y(n) + h*K3 )
　　 Runge-Kutta法是基于泰勒展开方法，因而需要所求解具有较好的光滑性。
　　 属性：差分方法
　　 《数值分析简明教程》-2 Editon -高等教育出版社- page 105 算法流程图
　　 
　　 代码维护：2005.6.14 DragonLord
　　**/
　　#include<iostream.h>
　　#include<stdio.h>
　　#include<math.h>
　　/*
　　举例方程：
　　 
　　y'= y - 2*x / y ( 0<x<1 )
　　 y(0) = 1
　　*/
　　double f(double x,double y)
　　{
　　double re;
　　if(x==0)re=1;
　　else re=y-2*x/y;
　　return re;
　　}
　　int main()
　　{
　　double x0,x1,y0,y1,h,k1,k2,k3,k4,y;
　　int N;
　　while(cin>>x0>>y0>>h>>N)
　　{
　　 int n=0;
　　 for(;n<N;n++)
　　 {
　　 x1=x0+h;
　　 y=sqrt(1+2*x1);
　　 k1=f(x0,y0);
　　 k2=f(x0+h/2,y0+h*k1/2);
　　 k3=f(x0+h/2,y0+h*k2/2);
　　 k4=f(x1,y0+h*k3);
　　 y1=y0+h*(k1+2*k2+2*k3+k4)/6;
　　 
　　 printf("%.1f %.4f %.4f\n",x1,y1,y);
　　 x0=x1;
　　 y0=y1;
　　 }
　　}
　　return 0;
　　}