龙格库塔法的编程
龙格库塔法的编程
#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));
x[0]=a;
y[0]=y0;
switch(style)
{
case 2:
for(i=0;i { x[i+1]=x[i]+h; k1=function(x[i],y[i]); k2=function(x[i]+h/2,y[i]+h*k1/2); y[i+1]=y[i]+h*k2; } break; case 3: for(i=0;i { x[i+1]=x[i]+h; k1=function(x[i],y[i]); k2=function(x[i]+h/2,y[i]+h*k1/2); k3=function(x[i]+h,y[i]-h*k1+2*h*k2); y[i+1]=y[i]+h*(k1+4*k2+k3)/6; } break; case 4: for(i=0;i { x[i+1]=x[i]+h; k1=function(x[i],y[i]); k2=function(x[i]+h/2,y[i]+h*k1/2); k3=function(x[i]+h/2,y[i]+h*k2/2); k4=function(x[i]+h,y[i]+h*k3); y[i+1]=y[i]+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(0 /* int main() { double x[6],y[6]; printf("用二阶龙格-库塔方法\n"); RungeKutta(1,0,1,5,x,y,2,function); for(int i=0;i<6;i++) printf("x[%d]=%f,y[%d]=%f\n",i,x[i],i,y[i]); printf("用三阶龙格-库塔方法\n"); RungeKutta(1,0,1,5,x,y,3,function); for(i=0;i<6;i++) printf("x[%d]=%f,y[%d]=%f\n",i,x[i],i,y[i]); printf("用四阶龙格-库塔方法\n"); RungeKutta(1,0,1,5,x,y,4,function); for(i=0;i<6;i++) printf("x[%d]=%f,y[%d]=%f\n",i,x[i],i,y[i]); return 1; 龙格库塔法的c++编程 2007-08-12 22:41:52| 分类：MatLab/Maple/Mat|字号订阅 from: https://www.wendangku.net/doc/f614804123.html,/forum/viewthread.php?tid=25801&highlight=%BF% E2%CB%FE 龙格库塔法的c++编程 CODE: #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)); x[0]=a; y[0]=y0; switch(style) { case 2: for(i=0;i { x[i+1]=x[i]+h; k1=function(x[i],y[i]); k2=function(x[i]+h/2,y[i]+h*k1/2); y[i+1]=y[i]+h*k2; } break; case 3: for(i=0;i { x[i+1]=x[i]+h; k1=function(x[i],y[i]); k2=function(x[i]+h/2,y[i]+h*k1/2); k3=function(x[i]+h,y[i]-h*k1+2*h*k2); y[i+1]=y[i]+h*(k1+4*k2+k3)/6; } break; case 4: for(i=0;i { x[i+1]=x[i]+h; k1=function(x[i],y[i]); k2=function(x[i]+h/2,y[i]+h*k1/2); k3=function(x[i]+h/2,y[i]+h*k2/2); k4=function(x[i]+h,y[i]+h*k3); y[i+1]=y[i]+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; } //例子求???=<<*-='1)10(/20y x y x y y /* int main() { double x[6],y[6]; printf("用二阶龙格-库塔方法\n"); RungeKutta(1,0,1,5,x,y,2,function); for(int i=0;i<6;i++) printf("x[%d]=%f,y[%d]=%f\n",i,x[i],i,y[i]); printf("用三阶龙格-库塔方法\n"); RungeKutta(1,0,1,5,x,y,3,function); for(i=0;i<6;i++) printf("x[%d]=%f,y[%d]=%f\n",i,x[i],i,y[i]); printf("用四阶龙格-库塔方法\n"); RungeKutta(1,0,1,5,x,y,4,function); for(i=0;i<6;i++) printf("x[%d]=%f,y[%d]=%f\n",i,x[i],i,y[i]); return 1; 龙格库塔求解微分方程数值解 工程中很多的地方用到龙格库塔求解微分方程的数值解， 龙格库塔是很重要的一种方法，尤其是四阶的，精确度相当的高。 此代码只是演示求一个微分方程： ???=∈*-='1]6.0,0[/20y x y x y y 的解，要求解其它的微分方程，可以自己定义借口函数， 退换程序里面的函数： float f(float , float)即可； CODE: /* /*编程者：刘艮平 /*完成日期：2004.5.29 /*e.g:y'=y-2x/y,x ∈[0,0.6] /* y(0)=1 /*使用经典四阶龙格-库塔算法进行高精度求解 /* y(i ＋1)＝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) */ #include #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 y[3]; //用于存放计算出的常微分方程数值解 int i=0; int j; //以下为函数的接口 printf("intput x0:"); scanf("%f",&x0); printf("input x1:"); scanf("%f",&x1); printf("input y[0]:"); scanf("%f",&y[0]); printf("input h:"); scanf("%f",&h); // 以下为核心程序 n=((x1-x0)/h); n=3; for(j=0;j { k1=h*f(x0,y[i]); //求K1 k2=h*f((x0+h/2),(y[i]+k1/2)); //求K2 k3=h*f((x0+h/2),(y[i]+k2/2)); //求K3 k4=h*f((x0+h),(y[i]+k3)); //求K4 y[i+1]=y[i]+((k1+2*k2+2*k3+k4)/6); //求yi+1 x0+=float(0.2); printf("y[%f]=%f\n",x0,y[i+1]); ++i; } } 龙格库塔法 一、基本原理： 龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于 此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。该算法是构建在数学支持的基础之上的。对于一阶精度的欧拉公式有： ???=*+=+),(1 11i i i i y x f k k h y y （1） 当用点i x 处的斜率近似值1k 与右端点1+i x 处的斜率2k 的算术平均值作为平均斜 率*k 的近似值，那么就会得到二阶精度的改进欧拉公式： ?????*++==+*+=+),() ,(2/)(12 1211k h y h x f k y x f k k k h y y i i i i i i （2） 依次类推，如果在区间[xi,xi+1]内多预估几个点上的斜率值K1、K2、……Km ，并用他们的加权平均数作为平均斜率K*的近似值，显然能构造出具有很高精度的高阶计算公式。经数学推导、求解，可以得出四阶龙格－库塔公式，也就是在工程中应用广泛的经典龙格－库塔算法： yi+1=yi+h*( K1+ 2*K2 +2*K3+ K4)/6 K1=f(xi,yi) K2=f(xi+h/2,yi+h*K1/2) K3=f(xi+h/2,yi+h*K2/2) K4=f(xi+h,yi+h*K3) 通常所说的龙格-库塔法是指四阶而言的，我们可以仿二阶、三阶的情形推导出常用的标准四阶龙格-库塔法公式 （1） 计算公式（1）的局部截断误差是 。 龙格-库塔法具有精度高，收敛，稳定（在一定条件下），计算过程中可以改变步长，不需要计算高阶导数等优点，但仍需计算 在一些点上的值，如四阶龙格-库塔法每计算一步需要计算四次 的值，这给实际计算带来一定的复杂性，因此，多用来计算“表头”。 二、小程序 #include #include #define f(x,y) (-1*(x)*(y)*(y)) void main(void) { double a,b,x0,y0,k1,k2,k3,k4,h; int n,i; printf("input a,b,x0,y0,n:"); scanf("%lf%lf%lf%lf%d",&a,&b,&x0,&y0,&n); printf("x0\ty0\tk1\tk2\tk3\tk4\n"); for(h=(b-a)/n,i=0;i!=n;i++) { k1=f(x0,y0); k2=f(x0+h/2,y0+k1*h/2); k3=f(x0+h/2,y0+k2*h/2); k4=f(x0+h,y0+h*k3); printf("%lf\t%lf\t",x0,y0); printf("%lf\t%lf\t",k1,k2); printf("%lf\t%lf\n",k3,k4); y0+=h*(k1+2*k2+2*k3+k4)/6; x0+=h; } printf("xn=%lf\tyn=%lf\n",x0,y0); } 运行结果： input a,b,x0,y0,n:0 5 0 2 20 x0 y0 k1 k2 k3 k4 0.000000 2.000000 -0.000000 -0.500000 -0.469238 -0.886131 0.250000 1.882308 -0.885771 -1.176945 -1.129082 -1.280060 0.500000 1.599896 -1.279834 -1.295851 -1.292250 -1.222728 0.750000 1.279948 -1.228700 -1.110102 -1.139515 -0.990162 1.000000 1.000027 -1.000054 -0.861368 -0.895837 -0.752852 1.250000 0.780556 -0.761584 -0.645858 -0.673410 -0.562189 1.500000 0.615459 -0.568185 -0.481668 -0.500993 -0.420537 1.750000 0.492374 -0.424257 -0.361915 -0.374868 -0.317855 2.000000 0.400054 -0.320087 -0.275466 -0.284067 -0.243598 2.250000 0.329940 -0.244935 -0.212786 -0.218538 -0.189482 2.500000 0.275895 -0.190295 -0.166841 -0.170744 -0.149563 2.750000 0.233602 -0.150068 -0.132704 -0.135399 -0.119703 3.000000 0.200020 -0.120024 -0.106973 -0.108868 -0.097048 3.250000 0.172989 -0.097256 -0.087300 -0.088657 -0.079618 3.500000 0.150956 -0.079757 -0.072054 -0.073042 -0.066030 3.750000 0.132790 -0.066124 -0.060087 -0.060818 -0.055305 4.000000 0.117655 -0.055371 -0.050580 -0.051129 -0.046743 4.250000 0.104924 -0.046789 -0.042945 -0.043363 -0.039833 4.500000 0.094123 -0.039866 -0.036750 -0.037072 -0.034202 4.750000 0.084885 -0.034226 -0.031675 -0.031926 -0.029571 xn=5.000000 yn=0.076927 我这还有一个程序，从别的论坛上看到的，你看看对你有没有帮助！ function change_step_RK(fun); % 4阶变步长龙格-库塔法解dy=fun(x,y)型方程% Example: % fun=inline('cos(x)+sin(y)'); % change_step_RK(fun); AbsT ol=1e-4; h=0.01; p=4; % p是对应的阶数 x0=0; % x0是起点 xN=1; % xN是终点 y0=0; % 初值 x=x0; y=y0; n=1; p21=2^p-1; while x(end) [x1,y1]=RK_f(fun,x(n),y(n),h); [x2,y2]=RK_f(fun,x(n),y(n),h/2); if abs(y1-y2)/p21 x2=x1; y2=y1; h=2*h; [x1,y1]=RK_f(fun,x(n),y(n),h); end else while abs(y1-y2)/p21>AbsTol; x1=x2; y1=y2; h=h/2; [x2,y2]=RK_f(fun,x(n),y(n),h/2); end end [xa,ya]=RK_f(fun,h,x(n),y(n)); x(n+1)=xa; y(n+1)=ya; n=n+1; end plot(x,y,'k'); function [xa,ya]=RK_f(fun,h,x,y); % 4阶龙格-库塔法解下一点的值 k1=fun(x,y); k2=fun(x+h/2,y+h*k1/2); k3=fun(x+h/2,y+h*k2/2); k4=fun(x+h,y+h*k3); xa=x+h; ya=y+h*(k1+k2*2+2*k3+k4)/6; 拉格朗日 function y=lagrange(x0,y0,x) n=length(x0);m=length(x); for i=1:m z=x(i); s=0.0; for k=1:n p=1.0; for j=1:n if j~=k p=p*(z-x0(j))/(x0(k)-x0(j)); end end s=p*y0(k)+s; end y(i)=s; end SOR迭代法的Matlab程序 function [x]=SOR_iterative(A,b) % 用SOR迭代求解线性方程组,矩阵A是方阵x0=zeros(1,length(b)); % 赋初值 tol=10^(-2); % 给定误差界 N=1000; % 给定最大迭代次数 [n,n]=size(A); % 确定矩阵A的阶 w=1; % 给定松弛因子 k=1; % 迭代过程 while k<=N x(1)=(b(1)-A(1,2:n)*x0(2:n)')/A(1,1); for i=2:n x(i)=(1-w)*x0(i)+w*(b(i)-A(i,1:i-1)*x(1:i-1)'-A(i,i+1:n)*x0(i+1:n)')/A(i,i); end if max(abs(x-x0))<=tol fid = fopen('SOR_iter_result.txt', 'wt'); fprintf(fid,'\n********用SOR迭代求解线性方程组的输出结果********\n\n'); fprintf(fid,'迭代次数: %d次\n\n',k); fprintf(fid,'x的值\n\n'); fprintf(fid, '%12.8f \n', x); break; end k=k+1; x0=x; end if k==N+1 fid = fopen('SOR_iter_result.txt', 'wt'); fprintf(fid,'\n********用SOR迭代求解线性方程组的输出结果********\n\n'); fprintf(fid,'迭代次数: %d次\n\n',k); fprintf(fid,'超过最大迭代次数，求解失败！'); fclose(fid); end Matlab中龙格-库塔(Runge-Kutta)方法原理及实现龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法。由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。该算法是构建在数学支持的基础之上的。龙格库塔方法的理论基础来源于泰勒公式和使用斜率近似表达微分，它在积分区间多预计算出几个点的斜率，然后进行加权平均，用做下一点的依据，从而构造出了精度更高的数值积分计算方法。如果预先求两个点的斜率就是二阶龙格库塔法，如果预先取四个点就是四阶龙格库塔法。一阶常微分方程可以写作：y'=f(x,y),使用差分概念。 (Yn+1-Yn)/h= f(Xn,Yn)推出（近似等于，极限为Yn'） Yn+1=Yn+h*f(Xn,Yn) 另外根据微分中值定理，存在0 Yn+1=Yn+h*f(Xn+th,Y(Xn+th)) 这里K＝f(Xn+th,Y(Xn+th))称为平均斜率，龙格库塔方法就是求得K的一种算法。 利用这样的原理，经过复杂的数学推导（过于繁琐省略），可以得出截断误差为O(h^5)的四阶龙格库塔公式： K1＝f(Xn,Yn); K2=f(Xn+h/2,Yn+(h/2)*K1); K3=f(Xn+h/2,Yn+(h/2)*K2); K4=f(Xn+h,Yn+h*K3); Yn+1=Yn+h*(K1+2K2+2K3+K4)*(1/6); 所以，为了更好更准确地把握时间关系，应自己在理解龙格库塔原理的基础上，编写定步长的龙格库塔函数，经过学习其原理，已经完成了一维的龙格库塔函数。 仔细思考之后，发现其实如果是需要解多个微分方程组，可以想象成多个微分方程并行进行求解，时间，步长都是共同的，首先把预定的初始值给每个微分方程的第一步，然后每走一步，对多个微分方程共同求解。想通之后发现，整个过程其实很直观，只是不停的逼近计算罢了。编写的定步长的龙格库塔计算函数： function [x,y]=runge_kutta1(ufunc,y0,h,a,b)%参数表顺序依次是微分方程组的函数名称，初始值向量，步长，时间起点，时间终点（参数形式参考了ode45函数） n=floor((b-a)/h);%求步数 x(1)=a;%时间起点 y(:,1)=y0;%赋初值，可以是向量，但是要注意维数 for ii=1:n x(ii+1)=x(ii)+h; k1=ufunc(x(ii),y(:,ii)); k2=ufunc(x(ii)+h/2,y(:,ii)+h*k1/2); k3=ufunc(x(ii)+h/2,y(:,ii)+h*k2/2); k4=ufunc(x(ii)+h,y(:,ii)+h*k3); y(:,ii+1)=y(:,ii)+h*(k1+2*k2+2*k3+k4)/6; %按照龙格库塔方法进行数值求解 end 调用的子函数以及其调用语句： function dy=test_fun(x,y) dy = zeros(3,1);%初始化列向量 dy(1) = y(2) * y(3); dy(2) = -y(1) + y(3); dy(3) = -0.51 * y(1) * y(2); 对该微分方程组用ode45和自编的龙格库塔函数进行比较，调用如下： [T,F] = ode45(@test_fun,[0 15],[1 1 3]); subplot(121) plot(T,F)%Matlab自带的ode45函数效果 title('ode45函数效果') [T1,F1]=runge_kutta1(@test_fun,[1 1 3],0.25,0,15);%测试时改变test_fun的函数维数，别忘记改变初始值的维数 subplot(122) plot(T1,F1)%自编的龙格库塔函数效果 title('自编的龙格库塔函数') matlab编的4阶龙格库塔法解微分方程的程序 2010-03-10 20:16 function varargout=saxplaxliu(varargin) clc,clear x0=0;xn=1.2;y0=1;h=0.1; [y,x]=lgkt4j(x0,xn,y0,h); n=length(x); fprintf(' i x(i) y(i)\n'); for i=1:n fprintf('%2d %4.4f %4.4f\n',i,x(i),y(i)); end function z=f(x,y) z=-2*x*y^2; function [y,x]=lgkt4j(x0,xn,y0,h) x=x0:h:xn; n=length(x); y1=x; y1(1)=y0; for i=1:n-1 K1=f(x(i),y1(i)); K2=f(x(i)+h/2,y1(i)+h/2*K1); K3=f(x(i)+h/2,y1(i)+h/2*K2); K4=f(x(i)+h,y1(i)+h*K3); y1(i+1)=y1(i)+h/6*(K1+2*K2+2*K3+K4); end y=y1; 结果： i x(i) y(i) 1 0.0000 1.0000 2 0.1000 0.9901 3 0.2000 0.9615 4 0.3000 0.9174 5 0.4000 0.8621 6 0.5000 0.8000 7 0.6000 0.7353 8 0.7000 0.6711 9 0.8000 0.6098 10 0.9000 0.5525 11 1.0000 0.5000 12 1.1000 0.4525 13 1.2000 0.4098 例题一 四阶龙格-库塔法的具体形式为： 1.1程序: /*e.g: y'=t-y,t∈[0,1] /*y(0)=0 /*使用经典四阶龙格-库塔算法进行高精度求解 /* y(i＋1)＝yi＋( K1+ 2*K2 +2*K3+ K4)/6 /* K1=h*f(ti,yi) /* K2=h*f(ti+h/2,yi+K1*h/2) /* K3=h*f(ti+h/2,yi+K2*h/2) /* K4=h*f(ti+h,yi+K3*h) */ #include float t; //范围下限 float h; //步长 float nn; //计算出的点的个数,即迭代步数 int n; //对nn取整 float k1,k2,k3,k4; float y[N]; //用于存放计算出的常微分方程数值解 float yy; //精确值 float error;//误差 int i=0,j; //以下为函数的接口 printf("input t0:"); scanf("%f",&t0); printf("input t:"); scanf("%f",&t); printf("input y[0]:"); scanf("%f",&y[0]); printf("input h:"); scanf("%f",&h); // 以下为核心程序 nn=(t-t0)/h; printf("nn=%f\n",nn); n=(int)nn; printf("n=%d\n",n); for(j=0;j<=n;j++) { yy=t0-1+exp(-t0); //解析解表达式 error=y[i]-yy; //误差计算 printf("y[%f]=%f yy=%f error=%f\n",t0,y[i],yy,error);//结果输出k1=h*f(t0,y[i]); //求K1 k2=h*f((t0+h/2),(y[i]+k1*h/2)); //求K2 k3=h*f((t0+h/2),(y[i]+k2*h/2)); //求K3 龙格库塔方法matlab实现~ function ff=rk(yy,x0,y0,h,a,b)%yy为y的导函数,x0,y0,为初值,h为步长,a,b为区间 c=(b-a)/h+1;i1=1; %c为迭代步数；i1为迭代步数累加值 y=y0;z=zeros(c,6); %z生成c行，5列的零矩阵存放结果； %每行存放c次迭代结果，每列分别存放k1~k4及y的结果 for x=a:h:b if i1<=c k1=feval(yy,x,y); k2=feval(yy,x+h/2,y+(h*k1)/2); k3=feval(yy,x+h/2,y+(h*k2)/2); k4=feval(yy,x+h,y+h*k3); y=y+(h/6)*(k1+2*k2+2*k3+k4); z(i1,1)=x;z(i1,2)=k1;z(i1,3)=k2;z(i1,4)=k3;z(i1,5)=k4;z(i1,6)=y; i1=i1+1; end end fprintf(‘结果矩阵，第一列为x(n),第二列~第五列为k1~k4，第六列为y（n+1）的结果') z %在命令框输入下列语句 %yy=inline('x+y'); %>> rk(yy,0,1,0.2,0,1) %将得到结果 %结果矩阵，第一列为x(n),第二列~第五列为k1~k4第六列为y（n+1）的结果 %z = % 0 1.0000 1.2000 1.2200 1.4440 1.2428 % 0.2000 1.4428 1.6871 1.7115 1.9851 1.5836 % 0.4000 1.9836 2.2820 2.3118 2.6460 2.0442 % 0.6000 2.6442 3.0086 3.0451 3.4532 2.6510 % 0.8000 3.4510 3.8961 3.9407 4.4392 3.4365 % 1.0000 4.4365 4.9802 5.0345 5.6434 4.4401 /** ***四阶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 for(;n 1. matlab 新建.m文件，编写龙格-库塔法求解函数 function [x,y]=runge_kutta1(ufunc,y0,h,a,b)%参数表顺序依次是微分方程组的函数名称，初始值向量，步长，时间起点，时间终点（参数形式参考了ode45函数） n=floor((b-a)/h); %求步数 x(1)=a;%时间起点 y(:,1)=y0;%赋初值，可以是向量，但是要注意维数 for ii=1:n x(ii+1)=x(ii)+h; k1=ufunc(x(ii),y(:,ii)); k2=ufunc(x(ii)+h/2,y(:,ii)+h*k1/2); k3=ufunc(x(ii)+h/2,y(:,ii)+h*k2/2); k4=ufunc(x(ii)+h,y(:,ii)+h*k3); y(:,ii+1)=y(:,ii)+h*(k1+2*k2+2*k3+k4)/6; %按照龙格库塔方法进行数值求解 end 2.另外再新建一个.,m文件，定义要求解的常微分方程函数 function dx=fun1(t,x) dx =zeros(2,1);%初始化列向量 dx(1) =0.08574*x(2)-1.8874*x(1)-20.17; dx(2) =1.8874*x(1)-0.08574*x(2)+115.16; 3,再新建一个.m文件，利用龙格-库塔方法求解常微分方程 [T1,F1]=runge_kutta1(@fun1,[46.30 1296 ],1,0,20); %求解步骤2定义的fun1常微分方程，@fun1是调用其函数指针，从0到20，间隔为1 subplot(122) plot(T1,F1)%自编的龙格库塔函数效果 title('自编的龙格库塔函数') grid 运行步骤3文件即可得到结果，F1为估测值 或者可以调用matlab自带函数ode45求解 方法如下： 龙格-库塔方法及其matlab实现 摘要：本文的目的数值求解微分方程精确解，通过龙格-库塔法，加以利用matlab为工具 达到求解目的。龙格-库塔(Runge-Kutta)方法是一种在工程上应用广泛的高精度单步算法，用于数值求解微分方程。MatLab软件是由美国Mathworks公司推出的用于数值计算和图形 处理的科学计算系统环境。MatLab是英文MATrix LABoratory（矩阵实验室）的缩写。在MratLab环境下，用户可以集成地进行程序设计、数值计算、图形绘制、输入输出、文件 管理等各项操作。 关键词：龙格-库塔 matlab 微分方程 1.前言 1.1：知识背景 龙格－库塔法（Runge-Kutta）是用于非线性常微分方程的解的重要的一类隐式或显式迭代法。这些技术由数学家卡尔·龙格和马丁·威尔海姆·库塔于1900年左右发明。通常所说的龙格库塔方法是相对四阶龙格库塔而言的，成为经典四阶龙格库塔法。该方法具有精度高，收敛，稳定，计算过程中可以改变步长不需要计算高阶导数等优点，但是仍需计算在 一些点上的值，比如四阶龙格-库塔法没计算一步需要计算四步，在实际运用中是有一定复杂性的。 Matlab是在20世纪七十年代后期的事：时任美国新墨西哥大学计算机科学系主任的Cleve Moler教授出于减轻学生编程负担的动机，为学生设计了一组调用LINPACK和EISPACK库 程序的“通俗易用”的接口，此即用FORTRAN编写的萌芽状态的MATLAB。 经几年的校际流传，在Little的推动下，由Little、Moler、Steve Bangert合作，于1984年成立了MathWorks公司，并把MATLAB正式推向市场。从这时起，MATLAB的内核 采用C语言编写，而且除原有的数值计算能力外，还新增了数据图视功能。 MATLAB以商品形式出现后，仅短短几年，就以其良好的开放性和运行的可靠性， 使原先控制领域里的封闭式软件包（如英国的UMIST，瑞典的LUND和SIMNON，德国的KEDDC）纷纷淘汰，而改以MATLAB为平台加以重建。在时间进入20世纪九十年代的时候，MATLAB已经成为国际控制界公认的标准计算软件。 到九十年代初期，在国际上30几个数学类科技应用软件中，MATLAB在数值计算方面独占 鳌头，而Mathematica和Maple则分居符号计算软件的前两名。Mathcad因其提供计算、 图形、文字处理的统一环境而深受中学生欢迎。 1.2研究的意义 精确求解数值微分方程，对龙格库塔的深入了解与正确运用，主要是在已知方程导数和初 值信息，利用计算机仿真时应用，省去求解微分方程的复杂过程。利用matlab强大的数值计算功能，省去认为计算的过程，达到快速精确求解数值微分方程。在实际生活中可以利 用龙格库塔方法和matlab的完美配合解决问题。 1.3研究的方法 对实例的研究对比，实现精度的要求，龙格库塔是并不是一个固定的公式，所以只是对典 型进行分析 2012-2013（1）专业课程实践论文二阶Runge-Kutta方法 董文峰，0818180123，R数学08-1班 由改进的Euler 方法得到： ()) ,(),(21121 211 K y h x hf K y x hf K K K y y i i i i i i ++==++=?????+ 凡满足条件式有一簇形如上式的计算格式，这些格式统称为二阶龙格—库塔格式。因此改进的欧拉格式是众多的二阶龙格—库塔法中的一种特殊格式。 若取1,0,2121212=== =c c b a ，就是另一种形式的二阶龙格 - 库塔公式。 ??????? ++==+=+)2,2() ,(12121K h y h x f K y x f K hK y y n n n n n n （1） 此计算公式称为变形的二阶龙格—库塔法。 二级龙格-库塔方法是显式单步式,每前进一步需要计算两个函数值。由上面的讨论可知,适当选择四个参数y0,a,b,n ，可使每步计算两次函数值的二阶龙格-库塔方法达到二阶精度。下面以式子（1）为依据利用VC++6.0编译程序进行问题的求解。 #include 龙格－库塔法（Runge-Kutta） 数值分析中，龙格－库塔法（Runge-Kutta）是用于模拟常微分方程的解的重要的一类隐式或显式迭代法。这些技术由数学家卡尔·龙格和马丁·威尔海姆·库塔于1900年左右发明。 经典四阶龙格库塔法 RK4””或者就是龙格库塔法的家族中的一个成员如此常用，以至于经常被称为“RK4 “龙格库塔法”。 令初值问题表述如下。 则，对于该问题的RK4由如下方程给出： 其中 这样，下一个值(yn+1)由现在的值(yn)加上时间间隔(h)和一个估算的斜率的乘积决定。该斜率是以下斜率的加权平均： k1是时间段开始时的斜率； k2是时间段中点的斜率，通过欧拉法采用斜率k1来决定y在点tn+h/2的值；k3也是中点的斜率，但是这次采用斜率k2决定y值； k4是时间段终点的斜率，其y值用k3决定。 当四个斜率取平均时，中点的斜率有更大的权值： RK4法是四阶方法，也就是说每步的误差是h5阶，而总积累误差为h4阶。 注意上述公式对于标量或者向量函数(y可以是向量)都适用。 显式龙格库塔法 显示龙格－库塔法是上述RK4法的一个推广。它由下式给出 其中 (注意：上述方程在不同著述中由不同但却等价的定义)。 要给定一个特定的方法，必须提供整数s(阶段数)，以及系数aij(对于1≤j 陕西科技大学 数值计算课程设计任务书 理学院信息与计算科学/应用数学专业信息08/数学08 班级学生： 题目：典型数值算法的C++语言程序设计 课程设计从2010 年 5 月17日起到2010 年 6 月18 日 1、课程设计的内容和要求（包括原始数据、技术要求、工作要求等）： 每人需作10个算法的程序、必做6题、自选4题。 对每个算法要求用C++语言进行编程。 必选题： 1、经典四阶龙格库塔法解一阶微分方程组 2、高斯列主元法解线性方程组 3、牛顿法解非线性方程组 4、龙贝格求积分算法 5、三次样条插值算法（压紧样条）用C++语言进行编程计算 依据计算结果，用Matlab画图并观察三次样条插值效果。 6、M次多项式曲线拟合，据计算结果，用Matlab画图并观察拟合效果。 自选题：自选4道其他数值算法题目.每道题目重选次数不得超过5次. 2、对课程设计成果的要求〔包括图表、实物等硬件要求〕： 1）提交课程设计报告 按照算法要求，用C++语言设计和开发应用程序，提交由算法说明；程序设计说明；系统技术文档(包括系统各模块主要流程图，软件测试方案与测试记录、软件调试和修改记录、测试结论、运行情况记录)，系统使用说明书，源程序代码为附录构成的课程设计报告。 2）课程设计报告版式要求 打印版面要求：A4纸，页边距：上2cm，下2cm，左2.5cm、右2cm；字体：正文宋体、小四号；行距：固定值20；页眉1.5cm ，页脚1.75cm；页码位于页脚居中打印；奇数页页眉“数值计算课程设计”，偶数页页眉“算法名称”，页眉宋体小5号；段落及层次要求：每节标题以四号黑体左起打印（段前段后各0.5行），节下为小节，以小四号黑体左起打印（段前段后各0.5行）。换行后以小四号宋体打印正文。节、小节分别以1、1.1、1.1.1依次标出，空一字符后接各部分的标题。 当论文结构复杂，小节以下的标题，左起顶格书写，编号依次用（1）、（2）……或1）、2）……顺序表示。字体为小四号宋体。 对条文内容采用分行并叙时，其编号用（a）、（b）……或a）、b）……顺序表示，如果编号及其后内容新起一个段落，则编号前空两个中文字符。3）设计报告装订顺序与规范 封面 数值计算课程设计任务书 目录 数值计算设计课程设计报告正文 设计体会及今后的改进意见 参考文献（资料） 左边缘装订 3 指导教师：日期： 教研室主任：日期： 龙格库塔法一、基本原理： 可以得出四阶龙格－库塔公式，也就是在工程中应用广泛的经典龙格－库塔算法，即： yi＋１＝yi＋h*( K1+ 2*K2 +2*K3+ K4)/6 K1=f(xi,yi) K2=f(xi+h/2,yi+h*K1/2) K3=f(xi+h/2,yi+h*K2/2) K4=f(xi+h,yi+h*K3) 通常所说的龙格-库塔法就是指四阶——龙格库塔法，我们可以仿二阶、三阶的情形推导出常用的标准四阶龙格-库塔法公式。 （1） 计算公式（1）的局部截断误差是。 龙格-库塔法具有精度高，收敛，稳定（在一定条件下），计算过程中可以改变步长，不需要计算高阶导数等优点，但仍需计算 在一些点上的值，如四阶龙格-库塔法每计算一步需要计算四次 的值，这给实际计算带来一定的复杂性，因此，多用来计算“表头”。 二、小程序 #include #define f(x,y) (-1*(x)*(y)*(y)) void main(void) { double a,b,x0,y0,k1,k2,k3,k4,h; int n,i; printf("input a,b,x0,y0,n:"); scanf("%lf%lf%lf%lf%d",&a,&b,&x0,&y0,&n); printf("x0\ty0\tk1\tk2\tk3\tk4\n"); for(h=(b-a)/n,i=0;i!=n;i++) { k1=f(x0,y0); k2=f(x0+h/2,y0+k1*h/2); k3=f(x0+h/2,y0+k2*h/2); k4=f(x0+h,y0+h*k3); printf("%lf\t%lf\t",x0,y0); printf("%lf\t%lf\t",k1,k2); printf("%lf\t%lf\n",k3,k4); y0+=h*(k1+2*k2+2*k3+k4)/6; x0+=h; } printf("xn=%lf\tyn=%lf\n",x0,y0); } 解微分方程的欧拉法，龙格-库塔法简单实例比较 欧拉方法（Euler method)用以对给定初值的常微分方程(即初值问题)求解分为前EULER 法、后退EULER 法、改进的EULER 法。 缺点： 欧拉法简单地取切线的端点作为下一步的起点进行计算，当步数增多时，误差会因积累而越来越大。因此欧拉格式一般不用于实际计算。 改进欧拉格式（向前欧拉公式）： 为提高精度，需要在欧拉格式的基础上进行改进。采用区间两端的斜率的平均值作为直线方程的斜率。改进欧拉法的精度为二阶。 算法： 微分方程的本质特征是方程中含有导数项，数值解法的第一步就是设法消除其导数值。对于常微分方程： (,)dy f x y dx = [,]x a b ∈ 0()y a y = 可以将区间[,]a b 分成n 段，那么方程在第i x 点有'()(,())i i i y x f x y x =，再用向前差商近似代替导数则为： ((1)())(,())i i i i y x y x f x y x h +-= 在这里，h 是步长，即相邻两个结点间的距离。因此可以根据i x 点和i y 的数值计算出1i y +来： 1(,)i i i i y y h f x y +=+?0,1,2,i L = 这就是向前欧拉公式。 改进的欧拉公式： 将向前欧拉公式中的导数(,)i i f x y 改为微元两端导数的平均，即上式便是梯形的欧拉公式。 可见，上式是隐式格式，需要迭代求解。为了便于求解，使用改进的欧拉公式： 数值分析中，龙格－库塔法（Runge-Kutta ）是用于模拟常微分方程的解的重要的一类隐式或显式迭代法。实际上，龙格-库塔法是欧拉方法的一种推广，向前欧拉公式将导数项简单取为(,)n n f x y ,而改进的欧拉公式将导数项取为两端导数的平均。 龙格-库塔方法的基本思想： 在区间1[,]n n x x +内多取几个点，将他们的斜率加权平均，作为导数的近似。 龙格库塔法的家族中的一个成员如此常用，以至于经常被称为“RK4”或者就是“龙格库塔法”。 令初值问题表述如下。 '(,)y f t y =00()y t y = 则，对于该问题的RK4由如下方程给出： 11234(22)6 n n h y y k k k k +=++++ 其中 1(,)n n k f t y = 21(,)22 n n h h k f t y k =++ 32(,)22 n n h h k f t y k =++ 43(,)n n k f t h y hk =++ 四阶龙格库塔法解微分方程 一、四阶龙格库塔法解一阶微分方程 ①一阶微分方程： cos y t ，初始值y(0)=0，求解区间[0 10]。 MATLAB 程序： %%%%%%%%%%% 四阶龙哥库塔法解一阶微分方程 %%%%%%%%%%% y'=cost %%%%%%%%%%% y(0)=0, 0≤t ≤10，h=0.01 %%%%%%%%%%% y=sint h=0.01; hf=10; t=0:h:hf; y=zeros(1,length(t)); y(1)=0; F=@(t,y)(cos(t)); for i=1:(length(t)-1) k1=F(t(i),y(i)); k2=F(t(i)+h/2,y(i)+k1*h/2); k3=F(t(i)+h/2,y(i)+k2*h/2); k4=F(t(i)+h,y(i)+k3*h); y(i+1)=y(i)+1/6*(k1+2*k2+2*k3+k4)*h; end subplot(211) plot(t,y,'-') xlabel('t'); ylabel('y'); title('Approximation'); span=[0,10]; p=y(1); [t1,y1]=ode45(F,span,p); subplot(212) plot(t1,y1) xlabel('t'); ylabel('y'); title('Exact'); 图1 ②一阶微分方程：()22*/x t x x t =- ，初始值x(1)=2，求解区间[1 3]。 MATLAB 程序： %%%%%%%%%%% 四阶龙哥库塔法解微分方程 %%%%%%%%%%% x'(t)=(t*x-x^2)/t^2 %%%%%%%%%%% x(1)=2, 1≤t ≤3, h=1/128 %%%%%%%%%%% 精确解：x(t)=t/(0.5+lnt) h=1/128; %%%%% 步长 tf=3; t=1:h:tf; x=zeros(1,length(t)); x(1)=2; %%%%% 初始值 F_tx=@(t,x)(t.*x-x.^2)./t.^2; for i=1:(length(t)-1) k_1=F_tx(t(i),x(i)); k_2=F_tx(t(i)+0.5*h,x(i)+0.5*h*k_1); k_3=F_tx((t(i)+0.5*h),(x(i)+0.5*h*k_2)); #include 函数功能 ode是专门用于解微分方程的功能函数，他有ode23,ode45,ode23s等等，采用的是Runge-Kutta算法。ode45表示采用四阶，五阶runge-kutta单步算法,截断误差为(Δx)3。解决的是Nonstiff(非刚性)的常微分方程.是解决数值解问题的首选方法,若长时间没结果,应该就是刚性的,换用ode23来解. 使用方法 [T,Y]=ode45(odefun,tspan,y0) odefun是函数句柄,可以是函数文件名,匿名函数句柄或内联函数名 tspan是区间[t0tf]或者一系列散点[t0,t1,...,tf] y0是初始值向量 T返回列向量的时间点 Y返回对应T的求解列向量 [T,Y]=ode45(odefun,tspan,y0,options) options是求解参数设置,可以用odeset在计算前设定误差,输出参数,事件等 [T,Y,TE,YE,IE]=ode45(odefun,tspan,y0,options) 在设置了事件参数后的对应输出 TE事件发生时间 YE事件解决时间 IE事件消失时间 sol=ode45(odefun,[t0tf],y0...) sol结构体输出结果 应用举例 1求解一阶常微分方程 程序: 一阶常微分方程 odefun=@(t,y)(y+3*t)/t^2;%定义函数 tspan=[14];%求解区间 y0=-2;%初值 [t,y]=ode45(odefun,tspan,y0); plot(t,y)%作图 title('t^2y''=y+3t,y(1)=-2,1 龙格库塔法的编程 #includematlab编的4阶龙格库塔法解微分方程的程序
四阶龙格库塔法的编程(赵)
龙格库塔方法matlab实现
四阶龙格库塔法原理C代码
龙格-库塔法MATLAB
龙格库塔方法及其matlab实现
二阶龙格库塔方法
龙格-库塔法(Runge-Kutta)matlab代码及含义
1、经典四阶龙格库塔法解一阶微分方程组
龙格库塔法-原理及程序实现
欧拉法与龙格库塔法比较分析
四阶龙格库塔法解微分方程
用龙格库塔法解微分方程组
MATLAB中龙格 库塔(RUNGE KUTTA)方法原理及实现
龙格库塔法的编程