《用matlab求解常微分方程》由会员分享，可在线阅读，更多相关《用matlab求解常微分方程（5页珍藏版）》请在金锄头文库上搜索。
1、实验六实验六用用 matlabmatlab 求解常微分方程求解常微分方程1 1微分方程的概念微分方程的概念未知的函数以及它的某些阶的导数连同自变量都由一已知方程联系在一起的方程称为微分方程。如果未知函数是一元函数，称为常微分方程。常微分方程的一般形式为F(t, y, y, y, y(n) 0如果未知函数是多元函数， 成为偏微分方程。 联系一些未知函数的一组微分方程组称为微分方程组。 微分方程中出现的未知函数的导数的最高阶解数称为微分方程的阶。 假设方程中未知函数及其各阶导数都是一次的，称为线性常微分方程，一般表示为y(n) a1(t)y(n1) an1(t)yan(t)y b(t)假设上式中的
2、系数ai(t),i 1,2,n均与t无关，称之为常系数。2 2常微分方程的解析解常微分方程的解析解dy y 1dt有些微分方程可直接通过积分求解.例如,一解常系数常微分方程可化为dy dtty ce 1.其中c为任意常数.有些常微分方程可用一些y 1,两边积分可得通解为技巧,如别离变量法,积分因子法,常数变异法,降阶法等可化为可积分的方程而求得解析解.线性常微分方程的解满足叠加原理,从而他们的求解可归结为求一个特解和相应齐次微分方程的通解.一阶变系数线性微分方程总可用这一思路求得显式解。高阶线性常系数微分方程可用特征根法求得相应齐次微分方程的基本解，再用常数变异法求特解。一阶常微分方程与高阶微
3、分方程可以互化，已给一个n阶方程(n1)y y, y y, y y12n设，可将上式化为一阶方程组y(n) f (t, y, y, y(n1)反过来，在许多情况下，一阶微分方程组也可化为高阶方程。所以一阶微分方程组与高阶常微分方程的理论与方法在许多方面是相通的， 一阶常系数线性微分方程组也可用特征根法求解。3 3微分方程的数值解法微分方程的数值解法除常系数线性微分方程可用特征根法求解，少数特殊方程可用初等积分法求解外，大部分微分方程无限世界，应用中主要依靠数值解法。考虑一阶常微分方程初值问题y1 y2y y23y ynn1yn f (t, y1, y2, yn)y(t) f (t, y(t),
4、t0 t tfy(t0) y0其中y (y1, y2, ym), f ( f1, f2, fm), y0 (y10, y20, ym0).所谓数值解法，就t t1 tn tfy ,k 0,1,n称是寻求y(t)在一系列离散节点0上的近似值khk tk1tk为步长，通常取为常量h。最简单的数值解法是 Euler 法。Euler 法的思路极其简单：在节点出用差商近似代替导数y(tk) 这样导出计算公式称为Euler 格式y(tk1) y(tk)hyk1 yk hf (tk, yk),k 0,1,2,他能求解各种形式的微分方程。Euler 法也称折线法。Euler 方法只有一阶精度，改良方法有二阶R
5、unge-Kutta 法、四阶Runge-Kutta 法、五阶Runge-Kutta-Felhberg 法和先行多步法等，这些方法可用于解高阶常微分方程组初值问题。边值问题采用不同方法，如差分法、 有限元法等。数值算法的主要缺点是它缺乏物理理解。4 4解微分方程的解微分方程的 MATLABMATLAB命令命令MATLAB中主要用 dsolve 求符号解析解，ode45,ode23,ode15s 求数值解。s=dsolve(方程 1, 方程 2,初始条件 1， 初始条件 2 ,自变量)用字符串方程表示，自变量缺省值为 t。导数用 D 表示，2 阶导数用D2 表示，以此类推。S 返回解析解。在方程
6、组情形，s 为一个符号结构。tout,yout=ode45(yprime,t0,tf,y0) 采用变步长四阶 Runge-Kutta 法和五阶 Runge-Kutta-Felhberg法求数值解，yprime 是用以表示 f(t,y)的 M 文件名，t0 表示自变量的初始值，tf 表示自变量的终值，y0 表示初始向量值。输出向量 tout 表示节点(t0,t1, ,tn)T,输出矩阵 yout 表示数值解，每一列对应 y 的一个分量。假设无输出参数，则自动作出图形。ode45 是最常用的求解微分方程数值解的命令，对于刚性方程组不宜采用。ode23 与ode45 类似， 只是精度低一些。 ode
7、12s 用来求解刚性方程组， 是用格式同 ode45。 可以用 helpdsolve, help ode45 查阅有关这些命令的详细信息.例例 1 1 求以下微分方程的解析解1y ay b2y sin(2x) y, y(0) 0, y(0) 13f f g, g g f , f (0) 1, g(0) 1方程1求解的 MATLAB代码为：clear;s=dsolve(Dy=a*y+b)结果为s =-b/a+exp(a*t)*C1方程2求解的 MATLAB代码为：clear;s=dsolve(D2y=sin(2*x)-y,y(0)=0,Dy(0)=1,x)simplify(s) %以最简形式显示
8、 s结果为s =(-1/6*cos(3*x)-1/2*cos(x)*sin(x)+(-1/2*sin(x)+1/6*sin(3*x)*cos(x)+5/3*sin(x)ans =-2/3*sin(x)*cos(x)+5/3*sin(x)方程3求解的 MATLAB代码为：clear;s=dsolve(Df=f+g,Dg=g-f,f(0)=1,g(0)=1)simplify(s.f) %s 是一个结构simplify(s.g)结果为ans =exp(t)*cos(t)+exp(t)*sin(t)ans =-exp(t)*sin(t)+exp(t)*cos(t)例例 求解微分方程y y t 1, y
9、(0) 1,先求解析解，再求数值解，并进行比较。由clear;s=dsolve(Dy=-y+t+1,y(0)=1,t)simplify(s)ty t e可得解析解为。下面再求其数值解，先编写M 文件%M 函数function f=fun8(t,y)f=-y+t+1;再用命令clear; close; t=0:0.1:1;y=t+exp(-t); plot(t,y); %化解析解的图形hold on; %保留已经画好的图形，如果下面再画图，两个图形和并在一起t,y=ode45(fun8,0,1,1);plot(t,y,ro); %画数值解图形，用红色小圈画xlabel(t),ylabel(y)结
10、果见图1.41.351.31.251.21.151.11.051y00.20.4t0.60.81图 解析解与数值解由图可见，解析解和数值解吻合得很好。例例 3 3 求方程ml mgsin,(0) 0,(0) 0的数值解.不妨取l 1,g 9.8,(0) 15.则上面方程可化为 9.8sin,(0) 15,(0) 0先看看有没有解析解.运行 MATLAB代码clear;s=dsolve(D2y=9.8*sin(y),y(0)=15,Dy(0)=0,t)simplify(s)知原方程没有解析解.下面求数值解.令y1, y2可将原方程化为如下方程组y1 y2y2 9.8sin(y1)y (0) 15
11、, y (0) 021建立 M 文件如下%M 文件function f=fun9(t,y)f=y(2), 9.8*sin(y(1); %f 向量必须为一列向量运行 MATLAB代码clear; close;t,y=ode45(fun9,0,10,15,0);plot(t,y(:,1); %画随时间变化图,y(:2)则表示 的值xlabel(t),ylabel(y1)结果见图16.516y115.515012345t678910图 7.2 数值解由图可见，随时间t周期变化。习题习题 16-616-61求以下微分方程的解析解2求方程的解析解和数值解,并进行比较3分别用 ode45 和 ode15s 求解 Van-del-Pol方程(1 x2)y 2xy, y(0) 1, y(0) 3d2xdx(1 x2) x 021000dtdtx(0) 0,x0) 1的数值解,并进行比较.
