四阶RungeKutta方法.docx
《四阶RungeKutta方法.docx》由会员分享，可在线阅读，更多相关《四阶RungeKutta方法.docx（13页珍藏版）》请在冰豆网上搜索。
四阶RungeKutta方法
实验题目3四阶Runge-Kutta方法
摘要
一阶常微分方程初值问题
（6.1）
的数值解法是近似计算中很重要的部分。
常微分方程初值问题的数值解法是求方程（6.1）的解在点列
上的近似值
，这里
是
到
的步长，一般略去下标记为
。
常微分方程初值问题的数值解法一般分为两大类：
（1）单步法：
这类方法在计算
时，只用到
、
和
，即前一步的值。
因此，在有了初值以后就可以逐步往下计算。
典型方法如龙格–库塔
方法。
（2）多步法：
这类方法在计算
时，除用到
、
和
以外，还要用
，即前面
步的值。
典型方法如Adams方法。
经典的
方法是一个四阶的方法，它的计算公式是：
（6.2）
方法的优点是：
单步法、精度高，计算过程便于改变步长，缺点是计算量较大，每前进一步需要计算四次函数值
。
前言
利用四阶龙格-库塔方法求解微分方程的初值问题
程序设计流程
问题1
（1）TestRK4（'ode1',1,[0-1],5,inline（'-x-1'））
TestRK4（'ode1',1,[0-1],10,inline（'-x-1'））
TestRK4（'ode1',1,[0-1],20,inline（'-x-1'））
（2）TestRK4（'ode2',1,[01],5,inline（'1./（x+1）'））
TestRK4（'ode2',1,[01],10,inline（'1./（x+1）'））
TestRK4（'ode2',1,[01],20,inline（'1./（x+1）'））
问题2
（1）TestRK4（'ode3',3,[10],5,inline（'x.^2.*（exp（x）-x）'））
TestRK4（'ode3',3,[10],10,inline（'x.^2.*（exp（x）-x）'））
TestRK4（'ode3',3,[10],20,inline（'x.^2.*（exp（x）-x）'））
（2）TestRK4（'ode4',3,[1-2],5,inline（'2*x./（1-2*x）'））
TestRK4（'ode4',3,[1-2],10,inline（'2*x./（1-2*x）'））
TestRK4（'ode4',3,[1-2],20,inline（'2*x./（1-2*x）'））
问题3
（1）TestRK4（'ode5',1,[01/3],5,inline（'x.^2+1/3*exp（-20*x）'））
TestRK4（'ode5',1,[01/3],10,inline（'x.^2+1/3*exp（-20*x）'））
TestRK4（'ode5',1,[01/3],20,inline（'x.^2+1/3*exp（-20*x）'））
（2）TestRK4（'ode6',1,[01],5,inline（'exp（-20*x）+sin（x）'））
TestRK4（'ode6',1,[01],10,inline（'exp（-20*x）+sin（x）'））
TestRK4（'ode6',1,[01],20,inline（'exp（-20*x）+sin（x）'））
（3）TestRK4（'ode7',1,[00],5,inline（'exp（x）.*sin（x）'））
TestRK4（'ode7',1,[00],10,inline（'exp（x）.*sin（x）'））
TestRK4（'ode7',1,[00],20,inline（'exp（x）.*sin（x）'））
实验所用函数
function[x,y]=RK4ODE（fun,xEnd,ini,h）
%RK4ODE用四阶Runge-Kutta法解初值问题dy/dx=f（x,y）,y（x0）=y0，在x处y的值
%
%Synopsis:
[x,y]=RK4ODE（fun,xEnd）
%[x,y]=RK4ODE（fun,xEnd,ini）
%[x,y]=RK4ODE（fun,xEnd,ini,h）
%
%Input:
fun=（string）初值问题的函数
%xEnd=使用Euler法的截止点
%ini=（optional）初始条件[x0y0]，默认为[00]
%h=（optional）步长，默认为0.05
%
%Output:
y=初值问题在x处y的近似值
ifnargin<3
ini=[00];%若未给初始条件，将初始条件设为[00]
end
ifnargin<4
h=0.05;%若未给出步长，将步长设为0.05
end
ini=ini（:
）;%将初始条件转为列向量，便于判断是否正确
[m,n]=size（ini）;
ifm~=2|n~=1
error（'初始值必须是一个含两个元素的向量[x0y0]'）;
end
x0=ini
（1）;%初始化x0
y0=ini
（2）;%初始化y0
x=（x0:
h:
xEnd）';%构建x向量
y=y0*ones（length（x）,1）;%初始化y向量
forj=2:
length（x）
k1=h*feval（fun,x（j-1）,y（j-1））;%三阶Runge-Kutta法的递推公式：
y（n+1）=y（n）+（k1+2*k2+2*k3+k4）/6
k2=h*feval（fun,x（j-1）+h/2,y（j-1）+k1/2）;%k1=h*f（x（n）,y（n））
k3=h*feval（fun,x（j-1）+h/2,y（j-1）+k2/2）;%k2=h*f（x（n）+h/2,y（n）+k1/2）
k4=h*feval（fun,x（j-1）+h,y（j-1）+k3）;%k3=h*f（x（n）+h/2,y（n）+k2/2）
y（j）=y（j-1）+（k1+2*k2+2*k3+k4）/6;%k4=h*f（x（n）+h,y（n）+k3）
end
functionTestRK4（fun,xEnd,ini,N,result）
h=（xEnd-ini
（1））/N;
[x,y]=RK4ODE（fun,xEnd,ini,h）;
y0=feval（result,x）;
plot（x,y,'ro',x,y0,'b-'）;
functiondydx=ode1（x,y）
dydx=x+y;
functiondydx=ode2（x,y）
dydx=-y.^2;
functiondydx=ode3（x,y）
dydx=2*y./x+x.^2.*exp（x）;
functiondydx=ode4（x,y）
dydx=（y.^2+y）./x;
functiondydx=ode5（x,y）
dydx=-20*（y-x.^2）+2*x;
functiondydx=ode6（x,y）
dydx=-20*y+20*sin（x）+cos（x）;
functiondydx=ode7（x,y）
dydx=-20*（y-exp（x）.*sin（x））+exp（x）.*（sin（x）+cos（x））;
思考题
1、实验一中
（1）数值解和解析解相同，
（2）数值解和解析解稍有不同，因为四阶Runge-Kutta方法是以小段的线性算法来近似获得微分方程的数值解，
（1）的准确解是1阶的，
（2）的准确解是无限阶的，因此对于
（1）数值解和解析解相同。
2、N越大，在固定长度的区间上取点就越多，步长就越小，因为四阶Runge-Kutta方法的误差是O（h^5），所以步长的减小有利于数值解的进一步精确。
3、N很小时，误差会变得非常大，这是因为准确解在计算区间内斜率变化非常快，如果N很小会使步长变大，造成非常大的误差
升级会员 