《解微分方程欧拉法,R-K法及其MATLAB实例.doc》由会员分享,可在线阅读,更多相关《解微分方程欧拉法,R-K法及其MATLAB实例.doc(4页珍藏版)》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、【精品文档】如有侵权,请联系网站删除,仅供学习与交流解微分方程欧拉法,R-K法及其MATLAB实例.精品文档.解微分方程的欧拉法,龙格-库塔法及其MATLAB简单实例欧拉方法(Euler method)用以对给定初值的常微分方程(即初值问题)求解分为前进EULER法、后退EULER法、改进的EULER法。缺点:欧拉法简单地取切线的端点作为下一步的起点进行计算,当步数增多时,误差会因积累而越来越大。因此欧拉格式一般不用于实际计算。改进欧拉格式:为提高精度,需要在欧拉格式的基础上进行改进。采用区间两端的斜率的平均值作为直线方程的斜率。改进欧拉法的精度为二阶。算法为:微分方程的本质特征是方程中含有导
2、数项,数值解法的第一步就是设法消除其导数值。对于常微分方程:xa,by(a) = y0可以将区间a,b分成n段,那么方程在第xi点有y(xi) = f(xi,y(xi),再用向前差商近似代替导数则为:在这里,h是步长,即相邻两个结点间的距离。因此可以根据xi点和yi点的数值计算出yi+1来:i=0,1,2,L这就是向前欧拉格式。改进的欧拉公式:将向前欧拉公式中的导数f(xi,yi)改为微元两端导数的平均,即上式便是梯形的欧拉公式。可见,上式是隐式格式,需要迭代求解。为了便于求解,使用改进的欧拉公式:数值分析中,龙格库塔法(Runge-Kutta)是用于模拟常微分方程的解的重要的一类隐式或显式迭
3、代法。实际上,龙格-库塔法是欧拉方法的一种推广,向前欧拉公式将导数项简单取为f(xn,yn),而改进的欧拉公式将导数项取为两端导数的平均。龙格-库塔方法的基本思想:在区间xn,xn+1内多取几个点,将他们的斜率加权平均,作为导数的近似。龙格库塔法的家族中的一个成员如此常用,以至于经常被称为“RK4”或者就是“龙格库塔法”。令初值问题表述如下。则,对于该问题的RK4由如下方程给出:其中这样,下一个值(yn+1)由现在的值(yn)加上时间间隔(h)和一个估算的斜率的乘积决定。该斜率是以下斜率的加权平均: k1是时间段开始时的斜率; k2是时间段中点的斜率,通过欧拉法采用斜率k1来决定y在点tn +
4、 h/2的值; k3也是中点的斜率,但是这次采用斜率k2决定y值; k4是时间段终点的斜率,其y值用k3决定。当四个斜率取平均时,中点的斜率有更大的权值:RK4法是四阶方法,也就是说每步的误差是h5阶,而总积累误差为h4阶。注意上述公式对于标量或者向量函数(y可以是向量)都适用。例子:下面给出了数值求解该微分方程的简单程序。其中y1,y2,y3,y4分别为向前欧拉公式,改进的欧拉公式,4级4阶龙格-库塔公式及精确解。h=0.1;x=0:h:1;y1=zeros(size(x);y1(1)=1;y2=zeros(size(x);y2(1)=1;y3=zeros(size(x);y3(1)=1;f
5、or i1=2:length(x) y1(i1)=y1(i1-1)+h*(y1(i1-1)-2*x(i1-1)/y1(i1-1); k1=y2(i1-1)-2*x(i1-1)/y2(i1-1); k2=y2(i1-1)+h*k1-2*x(i1)/(y2(i1-1)+h*k1); y2(i1)=y2(i1-1)+h*(k1+k2)/2; k1=y2(i1-1)-2*x(i1-1)/y2(i1-1); k2=y2(i1-1)+h*k1/2-2*(x(i1-1)+h/2)/(y2(i1-1)+h*k1/2); k3=y2(i1-1)+h*k2/2-2*(x(i1-1)+h/2)/(y2(i1-1)+h*k2/2); k4=y2(i1-1)+h*k3-2*(x(i1-1)+h)/(y2(i1-1)+h*k3); y3(i1)=y3(i1-1)+(k1+2*k2+2*k3+k4)*h/6; endy4=sqrt(1+2*x);%plot(x,y1,x,y2,x,y3,x,y4)%legend(y1,y2,y3,y4)plot(x,y4-y1,x,y4-y2,x,y4-y3)legend(y1,y2,y3)