常微分方程数值解法.doc

《常微分方程数值解法.doc》由会员分享，可在线阅读，更多相关《常微分方程数值解法.doc（4页珍藏版）》请在得力文库 - 分享文档赚钱的网站上搜索。

1、如有侵权，请联系网站删除，仅供学习与交流常微分方程数值解法【精品文档】第 4 页常微分方程数值解法【作用】微分方程建模是数学建模的重要方法，因为许多实际问题的数学描述将导致求解微分方程的定解问题。把形形色色的实际问题化成微分方程的定解问题，大体上可以按以下几步：1. 根据实际要求确定要研究的量(自变量、未知函数、必要的参数等)并确定坐标系。2. 找出这些量所满足的基本规律(物理的、几何的、化学的或生物学的等等)。3. 运用这些规律列出方程和定解条件。 基本模型1.发射卫星为什么用三级火箭2.人口模型3.战争模型4.放射性废料的处理通常需要求出方程的解来说明实际现象，并加以检验。如果能得到解析形

2、式的解固然是便于分析和应用的，但是我们知道，只有线性常系数微分方程，并且自由项是某些特殊类型的函数时，才可以得到这样的解，而绝大多数变系数方程、非线性方程都是所谓“解不出来”的于是对于用微分方程解决实际问题来说，数值解法就是一个十分重要的手段。1. 改进Euler 法：2. 龙格库塔（RungeKutta）方法：【源程序】1. 改进Euler 法：function x,y=eulerpro(fun,x0,x1,y0,n);%fun为函数，(x0，x1)为x区间，y0为初始值，n为子区间个数if nargin5,n=50;endh=(x1-x0)/n;x(1)=x0;y(1)=y0;for i=

3、1:nx(i+1)=x(i)+h;y1=y(i)+h*feval(fun,x(i),y(i);y2=y(i)+h*feval(fun,x(i+1),y1);y(i+1)=(y1+y2)/2;end调用 command窗口f=inline(-2*y+2*x2+2*x)x,y=eulerpro(f,0,0.5,1,10)求解函数y=2y+2+2x ，(0 x 0.5)， y(0) = 12. 龙格库塔（RungeKutta）方法：t,y=solver(F,tspan，y0)这里solver为ode45，ode23，ode113，输入参数 F 是用M文件定义的微分方程y= f (x, y)右端的函数

4、。tspan=t0,tfinal是求解区间，y0是初值。注：ode45和ode23变步长的，采用Runge-Kutta算法。ode45表示采用四阶-五阶Runge-Kutta算法，它用4阶方法提供候选解，5阶方法控制误差，是一种自适应步长（变步长）的常微分方程数值解法，其整体截断误差为(x)5。解决的是Nonstiff(非刚性)常微分方程。ode45是解决数值解问题的首选方法，若长时间没结果，应该就是刚性的，可换用ode23试试例如：odefun=(t,y) (y+3*t)/t2; %定义函数tspan=1 4; %求解区间y0=-2; %初值t,y=ode45(odefun,tspan,y0

5、);plot(t,y) %作图title(t2y=y+3t,y(1)=-2,1t4)legend(t2y=y+3t)xlabel(t)ylabel(y)% 精确解dsolve(t2*Dy=y+3*t,y(1)=-2) ans = (3*Ei(1,-1/t)+(-3*exp(-1)*Ei(1,-1)-2)/exp(-1)*exp(-1/t)求解高阶常微分方程需要求解的高阶常微分方程：求解的关键是将高阶转为一阶,odefun的书写.F(y,y,y.y(n-1),t)=0用变量替换,y1=y,y2=y.注意odefun方程定义为行向量程序:function Testode45tspan=3.9 4.

6、0; %求解区间y0=2 8; %初值t,x=ode45(odefun,tspan,y0);plot(t,x(:,1),-o,t,x(:,2),-*)legend(y1,y2)title(y =-t*y + et*y +3sin2t)xlabel(t)ylabel(y)endfunction y=odefun(t,x)y=zeros(2,1); % 列向量y(1)=x(2);y(2)=-t*x(1)+exp(t)*x(2)+3*sin(2*t);end【原理】改进Euler 法：式称为由Euler 公式和梯形公式得到的预测校正系统，也叫改进Euler 法。为了编程方便，将式子改写成龙格库塔（RungeKutta）方法：在区间 Xn, Xn+1 内多取几个点，将它们的斜率加权平均作为K ，就有可能构造出精度更高的计算公式。这就是龙格库塔方法的基本思想。