龙格库塔法RKF45Matlab实现
龙格库塔法RKF45的Matlab实现
2007-08-16 14:03:32| 分类：MatLab/Maple/Mat|字号订阅
4阶5级龙格库塔法用于解一阶微分方程（组），对于高阶微分方程，可以将其转换为一阶微分方程组求解。原程序由John.H.Mathews编写（数值方法matlab版），但只能解微分方程，不能解微分方程组。由LiuLiu@uestc修改，使之能够解微分方程组。该程序精度比matlab自带的ode45更高。
rkf45.m:
function [Rt Rx]=rkf45(f,tspan,ya,m,tol)
% Input:
% - f function column vector
% - tspan[a,b] left & right point of [a,b]
% - ya initial value column vector
% -m initial guess for number of steps
% -tol tolerance
% Output:
% - Rt solution: vector of abscissas
% - Rx solution: vector of ordinates
% Program by John.Mathews, improved by liuliu@uestc
if length(tspan)~=2
error('length of vector tspan must be 2.');
end
if ~isnumeric(tspan)
error('TSPAN should be a vector of integration steps.');
end
if ~isnumeric(ya)
error('Ya should be a vector of initial conditions.');
end
h = diff(tspan);
if any(sign(h(1))*h <= 0)
error('Entries of TSPAN are not in order.') ;
end
a=tspan(1);
b=tspan(2);
ya=ya(:);
a2 = 1/4; b2 = 1/4; a3 = 3/8; b3 = 3/32; c3 = 9/32; a4 = 12/13;
b4 = 1932/2197; c4 = -7200/2197; d4 = 7296/2197; a5 = 1;
b5 = 439/216; c5 = -8; d5 = 3680/513; e5 = -845/4104; a6 = 1/2;
b6 = -8/27; c6 = 2; d6 = -3544/2565; e6 = 1859/4104; f6 = -11/40;
r1 = 1/360; r3 = -128/4275; r4 = -2197/75240; r5 = 1/50;
r6 = 2/55; n1 = 25/216; n3 = 1408/2565; n4 = 2197/4104; n5 = -1/5;
big = 1e15;
h = (b-a)/m;
hmin = h/64;% 步长自适应范围下限
hmax = 64*h;% 步长自适应范围上限
max1 = 200;% 迭代次数上限
Y(1,:) = ya;
T(1) = a;
j = 1;
% tj = T(1);
br = b - 0.00001*abs(b);
while (T(j) if ((T(j)+h)>br), h = b - T(j); end %caculate values of k1...k6,y1...y6 tj = T(j); yj = Y(j,:); y1 = yj; k1 = h*feval(f,tj,y1); y2 = yj+b2*k1; if big k2 = h*feval(f,tj+a2*h,y2); y3 = yj+b3*k1+c3*k2; if big k3 = h*feval(f,tj+a3*h,y3); y4 = yj+b4*k1+c4*k2+d4*k3; if big k4 = h*feval(f,tj+a4*h,y4); y5 = yj+b5*k1+c5.*k2+d5*k3+e5*k4; if big y6 = yj+b6*k1+c6.*k2+d6*k3+e6*k4+f6*k5; if big k6 = h*feval(f,tj+a6*h,y6); err = abs(r1*k1+r3*k3+r4*k4+r5*k5+r6*k6); ynew = yj+n1*k1+n3*k3+n4*k4+n5*k5; % error and step size control if ( (err Y(j+1,:) = ynew; if ((tj+h)>br), T(j+1) = b; else T(j+1) = tj + h; end j = j+1; tj = T(j); end if (max(err)==0), s = 0; else s1 = 0.84*(tol.*h./err).^(0.25);% 最佳步长值 s=min(s1); end if ((s<0.75)&(h>2*hmin)), h = h/2; end if ((s>1.50)&(2*h if ( (big end % [Rt Rx]=[T' Y]; Rt=T'; Rx=Y; 使用方法： 首先编写方程（组）文件（注意与ode45不同，这儿方程组为1Xn数组： function dx= fun(t,x) dx=zeros(1,2); dx(1)=x(1)+x(2)*2+0*t; dx(2)=3*x(1)+x(2)*2+0*t; 然后使用： [Rt,Rx]=rkf45(@fun,[0,0.2],[6;4],100,1e-7) 龙格库塔法 龙格库塔法是常用于模拟常微分方程的解的重要的一类隐式或显式迭代法。这些技术由数学家C. Runge和M.W. Kutta于1900年左右发明。由于此算法精度高，采取措施对误差进行抑制，所以其实现原理也较复杂。 龙格库塔法是一种在工程上应用广泛的高精度单步算法，可以应用在物理、工程、控制、动力学中，如模糊控制、弹道分析以及分析光纤特性等，在系统仿真中得到广泛应用。 龙格库塔法源自于相应的泰勒级数方法，在每一插值节点用泰勒级数展开,其截断误差阶数也是，根据可省略更高阶的导数计算, 这种方法可构造任意阶数的龙格库塔法。其中4 阶龙格库塔法是最常用的一种方法。因为它相当精确、稳定、容易编程。在计算中一般不必使用高阶方法, 因为附加的计算误差可由增加精度来弥补。如果需要较高的精度, 可采取减小步长的方法即可。4 阶龙格库塔法的精度类似4 阶泰勒级数法的精度。 1、初值问题 对于一阶常微分方程的初值问题 根据常微分方程的理论可知，此初值问题的解在区间[a，b]上存在，且唯一。 2、离散化 取步长h=(b-a)/n，将区间[a , b]分成n个子区间： a=<=b 在其中任意两点的曲线段上，根据积分中值定理，一段光滑曲 线上至少有一点，它的斜率与整段曲线的平均斜率相同， 得=y’() (0<<1) 其中，= 可以将上式改写成y()=y()+h*K (2.1) 其中K为平均斜率，K=f() 公式（2.1）表明，如果能够确定平均斜率K，就可以根据（2.1）式得到y()的值。 欧拉法和龙格库塔法就是用不同方法确定不同精度的平均斜率K，从而求得y()的近似值。 3、Euler法 欧拉法虽然精度低，但它是最简单的一种显式单步法，也是龙 格库塔法的基础。 首先，令、为y() 及y()的近似值，并且令平均斜 率K=f()，即以点的斜率作为平均斜率K，便得到欧拉公式=+h* f() （3.1） 4、改进的欧拉法 此种方法是取、两点的斜率的平均值作为平均斜率K， 即K= ,其中、均为y()以及y()的近似值，就得到 改进后的欧拉公式（4.1） 第五章 常微分方程的数值解法 5.2 龙格-库塔法 一、教学目标及基本要求 通过对本节课的学习，使学生掌握常微分方程、常微分方程方程组的数值解法。 二、教学内容及学时分配 本节课主要介绍常微分方程的数值解法。具体内容如下： 讲授内容：龙格库塔方法、收敛性与稳定性 三、教学重点难点 1．教学重点：龙格库塔方法、收敛性与稳定性。 2. 教学难点：收敛性与稳定性。 四、教学中应注意的问题 多媒体课堂教学为主。适当提问，加深学生对概念的理解。 五、正文 龙格－库塔方法 引言 龙格－库塔方法的基本思想 '11()() ()()()(,())n n n n y x y x y y x y x hf y h ξξξ++-=?=+ 令*(,())K f y ξξ=称为区间1[,]n n x x +上的平均斜率，只要对*K 提供一种算法，即可推导出一种计算公式。 欧拉公式只是取点n x 的斜率作为区间1[,]n n x x +的平均斜率*K ，精度自然很低。 考察改进的欧拉公式 11211 ()22 n n y y h k k +=++1(,)n n k f x y =21(,)n n k f x h y hk =++ 它利用n x ，1n x +两点的斜率取算术平均，1n x +处斜率通过已知信息n y 用欧拉公式预报得到。 可以考虑设法在1[,]n n x x +多取几个点的斜率值，将它们加权平均作为区间 1[,]n n x x +的平均斜率*K 。这就是龙格－库塔方法的基本思想。 1、二阶龙格－库塔方法 考察1[,]n n x x +内一点,01n p n x x hp p +=+<≤，希望通过,n n p x x +两个点的斜率 12,K K 加权平均得到*K ，即令112((1))n n y y h K K λλ+=+-+ 取1(,)n n K f x y =，如何预报n p x +处斜率2K ？ 仿照改进的欧拉公式，先用欧拉公式预测()n p y x +的值n p y +： 1n p n y y phK +=+ 然后用n p y +计算2(,)n p n p K f x y ++=，从而得 112121[(1)](,)(,) n n n n n p n y y h K K K f x y K f x y phK λλ++=+-+==+ 适当选取,p λ，使上述公式具有较高得精度。假定()n n y y x =，分别将12 ,K K 泰勒展开： '1(,)()n n n K f x y y x == 212 ' '' 2 (,) (,)[(,)(,)(,)]()()()() n p n n n x n n n n y n n n n K f x y pkK f x y ph f x y f x y f x y O h y x phy x O h +=+=+++=++ 代入得 '2''31()()()()n n n n y y x hy x ph y x O h λ+=+++ 按泰勒展开2' '' 31()()()()2 n n n n h y y x hy x y x O h +=+++ 比较得，只要1 2 p λ= ，公式截断误差为3()O h 特别，当1,1/2p λ==，就是改进的欧拉公式， 改取1,1/2p λ==， 12n n y y hk +=+，1(,)n n k f x y =，21(,)22 n n h h k f x y k =++ 四阶龙格一库塔法 通常所说的龙格一库塔法是指四阶而言的．我们可以仿二阶、三阶的情形推导出常用的标准四阶龙格一库塔公式 (9.22) 公式(9.22)的局部截断误差的阶为． 龙格一库塔法具有精度高，收敛，稳定(在一定的条件下)，计算过程中可以改变步长，不需要计算高阶导数值等优点．但仍需计算在一些点上的值，如四阶龙格－库塔法每计算一步需要算四次 的值，这给实际计算带来一定的复杂性，因此，多用来计算“表头”．(即开始几点的近似值)．例3.用标准龙格一库塔法求初值问题 在处的解． 解因与．若应用标准龙格一库塔方法公式(9.22)计算，对于n=0时，则有 于是得 这个值与准确解在处的值已十分接近．再对n=1，2，3，4应用式(9.22)计算，具体计算结果如表3所示： 例3写出用四阶龙格――库塔法求解初值问题 的计算公式，取步长h＝0.2计算y(0.4)的近似值。至少保留四位小数。 解此处f(x，y)＝8－3y，四阶龙格――库塔法公式为 其中κ1＝f(x k，y k)；κ2＝f(x k+0.5h，y k+0.5hκ1)； κ3＝f(x k+0.5h，y k+0.5hκ2)；κ4＝f(x k+h，y k+hκ3) 本例计算公式为： 其中κ1＝8－3y k；κ2＝5.6－2.1y k； κ3＝6.32－2.37y k；κ4＝4.208－1.578y k ＝1.2016＋0.5494y k (k＝0，1，2，…) 当x0＝0，y0＝2， y(0.2)≈y1＝1.2016＋0.5494y0＝1.2016＋0.5494×2＝2.3004 y(0.4)≈y2＝1.2016＋0.5494y1＝1.2016＋0.5494×2.3004＝2.4654 常微分方程组的四阶Runge-Kutta方法1.问题： 1.1若用普通方法-----仅适用于两个方程组成的方程组 编程实现： 创建M 文件： function R = rk4(f,g,a,b,xa,ya,N) %UNTITLED2 Summary of this function goes here % Detailed explanation goes here %x'=f(t,x,y) y'=g(t,x,y) %N为迭代次数 %h为步长 %ya,xa为初值 f=@(t,x,y)(2*x-0.02*x*y); g=@(t,x,y)(0.0002*x*y-0.8*y); h=(b-a)/N; T=zeros(1,N+1); X=zeros(1,N+1); Y=zeros(1,N+1); T=a:h:b; X(1)=xa; Y(1)=ya; for j=1:N f1=feval(f,T(j),X(j),Y(j)); g1=feval(g,T(j),X(j),Y(j)); f2=feval(f,T(j)+h/2,X(j)+h/2*f1,Y(j)+g1/2); g2=feval(g,T(j)+h/2,X(j)+h/2*f1,Y(j)+h/2*g1); f3=feval(f,T(j)+h/2,X(j)+h/2*f2,Y(j)+h*g2/2); g3=feval(g,T(j)+h/2,X(j)+h/2*f2,Y(j)+h/2*g2); f4=feval(f,T(j)+h,X(j)+h*f3,Y(j)+h*g3); g4=feval(g,T(j)+h,X(j)+h*f3,Y(j)+h*g3); X(j+1)=X(j)+h*(f1+2*f2+2*f3+f4)/6; Y(j+1)=Y(j)+h*(g1+2*g2+2*g3+g4)/6; R=[T' X' Y']; end 情况一：对于x0=3000，y0=120 控制台中输入： >> rk4('f','g',0,10,3000,120,10) 运行结果： ans = 1.0e+003 * 0 3.0000 0.1200 0.0010 2.6637 0.0926 0.0020 3.7120 0.0774 0.0030 5.5033 0.0886 0.0040 4.9866 0.1193 0.0050 3.1930 0.1195 0.0060 2.7665 0.0951 0.0070 3.6543 0.0799 0.0080 5.2582 0.0884 0.0090 4.9942 0.1157 0.0100 3.3541 0.1185 数据: 实验五变步长的龙格库塔法 一、实验目的 1、了解变步长的龙哥库塔法的特点及具体实现过程。 2、运用vc6.0软件编写龙哥库塔法的算法程序。 3、采用编制的程序，对实际问题进行数值模拟计算。 4、将数值计算结果与精确解进行对比分析，讨论计算误差。 二、实验内容 1、根据变步长的龙哥库塔法法算法的特点，设计程序的流程 本实验采用经典的变步长的龙哥库塔法法格式（四阶变步长的龙哥库塔法法格式） 2、用vc6.0语言编写变步长的龙哥库塔法法算法程序 在编写程序时，充分考虑到程序的交互性和实用性。本程序可以实现计算机主动提示操作步骤和输出结果； 输入待计算的信息后，计算机自动生成图像，并且将得到的结果也显示在图像中，便于直观观察； 3、算法的源程序 #include "stdio.h" #include "math.h" double f(double,double); main() { double y0=1,k1,k2,k3,k4,x0=0,h=0.1,h1=0.05,q,tol=0.1,y00,y01; int i,n=10; k1=f(x0,y0); k2=f(x0+h/2,y0+h/2*k1); k3=f(x0+h/2,y0+h/2*k2); k4=f(x0+h/2,y0+h*k3); y0=y0+h/6*(k1+2*k2+2*k3+k4); y00=y0; printf("the y00 is: %lf\n",y00); /*步长为h时从x0出发求一步得y00*/ k1=f(x0,y0); k2=f(x0+h1/2,y0+h1/2*k1); k3=f(x0+h1/2,y0+h1/2*k2); k4=f(x0+h1/2,y0+h1*k3); y0=y0+h1/6*(k1+2*k2+2*k3+k4); x0=x0+h1; k1=f(x0,y0); k2=f(x0+h1/2,y0+h1/2*k1); k3=f(x0+h1/2,y0+h1/2*k2); 龙格库塔法一、基本原理： 可以得出四阶龙格－库塔公式，也就是在工程中应用广泛的经典龙格－库塔算法，即： 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); } 龙格库塔方法 从常微分方程数值解法的几何意义看，欧拉方法取一点n t 处的斜率() n n Q t f k ,1 =作为平均斜率，因此欧拉方法近似公式为 t k Q Q n n ?+=+11 向后的欧拉方法则采用点1+n t 处的斜率 ()112,++=n n Q t f k 作为平均斜率，即 t k Q Q n n ?+=+21 所以这两种方法也称作矩形法。改进的欧拉方法则取点n t 处和点1+n t 处斜率 1k 和2k 的平均值作为平均斜率，即 ()t k k Q Q n n ?++=+2112 1 因此改进的欧拉方法又称为梯形方法。可以预见，若去多点处斜率的加权平均值作为平均斜率，误差会更小，这就是龙格-库塔方法。 最常用的是四阶龙格库卡近似计算公式，即： ()t k k k k Q Q n n ?++++=+43211226 1 式中 () () 314221312 121,2,2,,tk Q t f k k t Q t f k k t Q t f k Q t f k n n n n n n n n ?+=??? ? ???+=???? ???+==+++ 使用标准四阶龙格——库塔法求解初值问题 ()(){}10210≤≤='=x xy y y 计算所用程序如下所示 #include "stdio.h" #include "conio.h" float func(float x,float y) { return(2*x*y); } float runge_kutta(float x0,float xn,float y0,int n) { float x,y,y1,y2,h,xh; float d1,d2,d3,d4; int i; x=x0; y=y0; h=(xn-x0)/n; for(i=1;i<=n;i++) { xh=x+h/2; d1=func(x,y); d2=func(xh,y+h*d1/2.0); d3=func(xh,y+h*d2/2.0); d4=func(xh,y+h*d3); y=y+h*(d1+2*d2+2*d3+d4)/6.0; x=x0+i*h; } return(y); } void main() { float x0,xn,y0,e; int n; printf("\ninput n:\n"); scanf("%d",&n); printf("input x0,xn:\n"); scanf("%f%f",&x0,&xn); printf("input y0:\n"); scanf("%f",&y0); e=runge_kutta(x0,xn,y0,n); 8.2 龙格-库塔方法 8.2.1 二阶龙格-库塔方法 常微分方程初值问题： 做在点的泰勒展开： 这里。取，就有 (8.11) 截断可得到近似值的计算公式，即欧拉公式： 若取，式（8.11）可写成： 或 （8.12） 截断可得到近似值的计算公式： 或 上式为二阶方法，一般优于一阶的欧拉公式（8.2），但是在计算时，需要计算 在点的值，因此，此法不可取。 龙格-库塔设想用在点和值的线性组合逼近式（8.12）的主体，即用 （8.13） 逼近 得到数值公式： （8.14） 或更一般地写成 对式（8.13）在点泰勒展开得到： 将上式与式（8.12）比较，知当满足 时有最好的逼近效果，此时式（8.13）－式（8.14）。这是4个未知数的3个方程，显然方程组有无数组解。 若取，则有二阶龙格-库塔公式，也称为改进欧拉公式： （8.15） 若取，则得另一种形式的二阶龙格-库塔公式，也称中点公式： （8.16） 从公式建立过程中可看到，二阶龙格-库塔公式的局部截断误差仍为，是二阶精度的计算公式。类似地，可建立高阶的龙格-库塔公式，同时可知四阶龙格-库塔公式的局部 截断误差为，是四阶精度的计算公式。 欧拉法是低精度的方法，适合于方程的解或其导数有间断的情况以及精度要求不高的情况，当解需要高精度时，必须用高阶的龙格-库塔等方法。 四阶龙格-库塔方法应用面较广，具有自动起步和便于改变步长的优点，但计算量比一般方法略大。为了保证方法的收敛性，有时需要步长取得较小，因此，不适于解病态方程。 8.2.2 四阶龙格-库塔公式 下面列出常用的三阶、四阶龙格-库塔计算公式。 三阶龙格-库塔公式 龙格库塔法的编程 #include 利用龙格-库塔法求解质点运动常微分方程 一、待解问题 讨论常微分方程的初值问题，边值问题的数值解法，最常用的基本方法就是龙格-库塔法使一些物理方程的计算简便，所得结果的精确提高。 下面是利用龙格-库塔法求解一个质点运动的常微分方程 以初速s km v /80=自地球表面（半径km R 6000=）竖直向上发射一质量为 m 的火箭，如图所示。若不计空气阻力，火箭所受引力 F 大小与它到地心距离的平方成反比，求火箭所能到达的最大高度H 。 火箭为我们分析的对象，对其做受力分析，火箭在任意位置 x 处仅受地球引力F 作用。由题意知，F 的大小与 x 2 成反比，设μ为比例系数，则有： （1） 当火箭处于地面时，即R x =时F = m g ，由式(1)可得2mgR =μ，于 是火箭在任意位置 x 处所受地球引力 F 的大小为 22 x mgR F = 由于火箭作直线运动，火箭的直线运动微分方程式为： 2x F μ =2 222x mgR dt x d m -= 分离变量积分 22 vx gR dx dv -= 二、物理机理 1.龙格-库塔法的基本思想及一般形式 设初值问题],[)(b a x y y ∈=，由微分中值定理可知，必存在],[1+∈n n x x ξ，使 设)(n n x y y =，并记))(,(*ξξy f K =，则 *1)(hK y x y n n +=+ 其中*K 称为)(x y 在][1,+n n x x 上的平均斜率，只要对平均斜率*K 提供一种算法，上式就给出了一种数值解公式，例如，用),(1n n y x f K =代替*K ，就得到欧拉公式，用),(112++=n n y x f K 代替*K ，就得到向后欧拉公式，如果用21,K K 的平均值来代替*K ，则可得到二阶精度的梯形公式。可以设想，如果在],[1+n n x x 上能多预测几个点的斜率值，用它们的加权平均值代替*K ，就有望的到具有较高精度的数值解公式，这就是龙格-库塔（Runge-Kutta ）法的基本思想。 龙格-库塔公式的一般形式： )) (,()()()()(1ξξξy hf x y y h x y x y n n n +='+=+dx dv v dt dx dx dv dt dv dt x d =?==2 22 2 x mgR dx dv mv - =?龙格库塔积分算法
5.2龙格库塔法
龙格库塔法例题
常微分方程组的四阶RungeKutta龙格库塔法matlab实现
实验五 变步长的龙哥库塔法
龙格库塔法-原理及程序实现
龙格库塔方法解题
龙格-库塔方法
龙格库塔法的编程
利用龙格库塔法求解质点运动方程