/** ***四阶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; }
|