《龙格库塔法的编程》由会员分享，可在线阅读，更多相关《龙格库塔法的编程（16页珍藏版）》请在金锄头文库上搜索。
1、龙格库塔法的编程#include #include /*n表示几等分，n+1表示他输出的个数*/ int RungeKutta(double y0,double a,double b,int n,double *x,double *y,int style,double (*function)(double，double) double h=(b-a)/n,k1,k2,k3,k4; int i; / x=(double*)malloc(n+1)*sizeof(double); / y=(double*)malloc(n+1)*sizeof(double); x0=a; y0=y0; switch
2、(style) case 2: for(i=0;in;i+) xi+1=xi+h; k1=function(xi,yi); k2=function(xi+h/2,yi+h*k1/2); yi+1=yi+h*k2; break; case 3: for(i=0;in;i+) xi+1=xi+h; k1=function(xi,yi); k2=function(xi+h/2,yi+h*k1/2); k3=function(xi+h,yi-h*k1+2*h*k2); yi+1=yi+h*(k1+4*k2+k3)/6; break; case 4: for(i=0;in;i+) xi+1=xi+h;
3、k1=function(xi,yi);k2=function(xi+h/2,yi+h*k1/2); k3=function(xi+h/2,yi+h*k2/2); k4=function(xi+h,yi+h*k3);yi+1=yi+h*(k1+2*k2+2*k3+k4)/6; break; default: return 0; return 1;double function(double x,double y)return y-2*x/y;/例子求y=y-2*x/y(0x1);y0=1;/* int main() double x6,y6;printf(用二阶龙格-库塔方法n);RungeKu
4、tta(1,0,1,5,x,y,2,function);for(int i=0;i6;i+)printf(x%d=%f,y%d=%fn,i,xi,i,yi);printf(用三阶龙格-库塔方法n);RungeKutta(1,0,1,5,x,y,3,function);for(i=0;i6;i+)printf(x%d=%f,y%d=%fn,i,xi,i,yi); printf(用四阶龙格-库塔方法n); RungeKutta(1,0,1,5,x,y,4,function); for(i=0;i6;i+) printf(x%d=%f,y%d=%fn,i,xi,i,yi); return 1;龙格库
5、塔法的c+编程2007-08-12 22:41:52|分类：MatLab/Maple/Mat|字号订阅from:http:/ RungeKutta(double y0,double a,double b,int n,double *x,double *y,int style,double (*function)(double,double)double h=(b-a)/n,k1,k2,k3,k4;int i;/ x=(double*)malloc(n+1)*sizeof(double);/ y=(double*)malloc(n+1)*sizeof(double);x0=a;y0=y0;swi
6、tch(style)case 2:for(i=0;in;i+)xi+1=xi+h;k1=function(xi,yi);k2=function(xi+h/2,yi+h*k1/2);yi+1=yi+h*k2;break;case 3:for(i=0;in;i+)xi+1=xi+h;k1=function(xi,yi);k2=function(xi+h/2,yi+h*k1/2);k3=function(xi+h,yi-h*k1+2*h*k2);yi+1=yi+h*(k1+4*k2+k3)/6;break;case 4:for(i=0;in;i+)xi+1=xi+h;k1=function(xi,y
7、i);k2=function(xi+h/2,yi+h*k1/2);k3=function(xi+h/2,yi+h*k2/2);k4=function(xi+h,yi+h*k3);yi+1=yi+h*(k1+2*k2+2*k3+k4)/6;break;default:return 0;return 1;double function(double x,double y)return y-2*x/y;/例子求/*int main()double x6,y6;printf(用二阶龙格-库塔方法n);RungeKutta(1,0,1,5,x,y,2,function);for(int i=0;i6;i
8、+)printf(x%d=%f,y%d=%fn,i,xi,i,yi);printf(用三阶龙格-库塔方法n);RungeKutta(1,0,1,5,x,y,3,function);for(i=0;i6;i+)printf(x%d=%f,y%d=%fn,i,xi,i,yi);printf(用四阶龙格-库塔方法n);RungeKutta(1,0,1,5,x,y,4,function);for(i=0;i6;i+)printf(x%d=%f,y%d=%fn,i,xi,i,yi);return 1;龙格库塔求解微分方程数值解工程中很多的地方用到龙格库塔求解微分方程的数值解，龙格库塔是很重要的一种方法，
9、尤其是四阶的，精确度相当的高。此代码只是演示求一个微分方程：的解，要求解其它的微分方程，可以自己定义借口函数，退换程序里面的函数：float f(float , float)即可；CODE:/*/*编程者：刘艮平/*完成日期：2004.5.29/*e.g:y=y-2x/y,x0,0.6/*y(0)=1/*使用经典四阶龙格-库塔算法进行高精度求解/* y(i1)yi( K1+ 2*K2 +2*K3+ K4)/6/* K1=h*f(xi,yi)/* K2=h*f(xi+h/2,yi+K1/2)/* K3=h*f(xi+h/2,yi+K2/2)/* K4=h*f(xi+h,yi+K3)*/#incl
10、ude #include float f(float den,float p0)/要求解的微分方程的右部的函数 e.g: y-2x/y float rus;/ den=w/(W0+sl);/ rus=k*A*f/Wp*sqrt(RTp)-(k-1)*. rus=p0-2*den/p0; return(rus);/float fden()/void main() float x0; /范围上限 int x1; /范围下限 float h;/步长 int n; /计算出的点的个数 float k1,k2,k3,k4; float y3; /用于存放计算出的常微分方程数值解 int i=0; int
11、 j;/以下为函数的接口 printf(intput x0:); scanf(%f,&x0); printf(input x1:); scanf(%f,&x1); printf(input y0:); scanf(%f,&y0); printf(input h:); scanf(%f,&h);/ 以下为核心程序 n=(x1-x0)/h); n=3; for(j=0;jn;j+) k1=h*f(x0,yi); /求K1 k2=h*f(x0+h/2),(yi+k1/2); /求K2 k3=h*f(x0+h/2),(yi+k2/2); /求K3 k4=h*f(x0+h),(yi+k3); /求K4
12、yi+1=yi+(k1+2*k2+2*k3+k4)/6); /求yi+1 x0+=float(0.2); printf(y%f=%fn,x0,yi+1); +i; 龙格库塔法一、基本原理：龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。该算法是构建在数学支持的基础之上的。对于一阶精度的欧拉公式有： （1）当用点处的斜率近似值与右端点处的斜率的算术平均值作为平均斜率的近似值，那么就会得到二阶精度的改进欧拉公式： （2）依次类推，如果在区间xi,xi+1内多预估几个点上的斜率值K1、K2、Km，并用他们的加权平均数作为平均斜率K*的近似值，显然能构造出具有很高精度的高阶计算公式。经数学推导、求解，可以得出四阶龙格库塔公式，也就
