下载文档原格式
matlab解常微分方程例题当涉及到使用MATLAB解常微分方程（ODE）的例题时，我们可以采用MATLAB中的ode45函数来进行求解。
ode45是一种常用的ODE求解器，它基于龙格-库塔方法。
下面我将以一个简单的例题来说明如何使用MATLAB解常微分方程。
假设我们要解决以下的常微分方程：dy/dt = -2y + 4t.初始条件为y(0) = 1。
首先，我们需要定义一个匿名函数来表示方程右侧的表达式，即-2y + 4t。
在MATLAB中，可以这样定义这个函数：f = @(t, y) -2y + 4t.接下来，我们需要定义时间范围和初始条件：tspan = [0 5] % 时间范围从0到5。
y0 = 1 % 初始条件y(0) = 1。
然后，我们可以使用ode45函数进行求解：[t, y] = ode45(f, tspan, y0)。
最后，我们可以绘制出解的图像：plot(t, y)。
xlabel('t')。
ylabel('y')。
title('Solution of dy/dt = -2y + 4t')。
这样，我们就得到了常微分方程的数值解，并用图像表示出来。
需要注意的是，这只是一个简单的例题，实际应用中可能会涉及更复杂的常微分方程。
但是使用MATLAB的ode45函数求解常微分方程的基本步骤是相似的，定义方程右侧的函数，设定时间范围和初始条件，然后使用ode45函数进行求解，并绘制出解的图像。
希望以上的解答能够满足你的需求。
如果你有更多关于MATLAB 解常微分方程的问题，欢迎继续提问。
MATLAB常微分⽅程的数值解法MATLAB常微分⽅程的数值解法⼀、实验⽬的科学技术中常常要求解常微分⽅程的定解问题，所谓数值解法就是求未知函数在⼀系列离散点处的近似值。
⼆、实验原理三、实验程序1. 尤拉公式程序四、实验内容选⼀可求解的常微分⽅程的定解问题，分别⽤以上1, 4两种⽅法求出未知函数在节点处的近似值，并对所求结果与分析解的（数值或图形）结果进⾏⽐较。
五、解答1. 程序求解初值问题取n=10源程序：euler23.m:function [A1,A2,B1,B2,C1,C2]=euler23(a,b,n,y0)%欧拉法解⼀阶常微分⽅程%初始条件y0h = (b-a)/n; %步长h%区域的左边界a%区域的右边界bx = a:h:b;m=length(x);%前向欧拉法y = y0;for i=2:my(i)=y(i-1)+h*oula(x(i-1),y(i-1));A1(i)=x(i);A2(i)=y(i);endplot(x,y,'r-');hold on;%改进欧拉法y = y0;for i=2:my(i)=y(i-1)+h/2*( oula(x(i-1),y(i-1))+oula(x(i),y(i-1))+h*(oula(x(i-1),x(i-1))));B1(i)=x(i);B2(i)=y(i);endplot(x,y,'m-');hold on;%欧拉两步公式y=y0;y(2)=y(1)+h*oula(x(1),y(1));for i=2:m-1y(i+1)=y(i-1)+2*h*oula(x(i),y(i));C1(i)=x(i);C2(i)=y(i);endplot(x,y,'b-');hold on;%精确解⽤作图xx = x;f = dsolve('Dy=-3*y+8*x-7','y(0)=1','x');%求出解析解y = subs(f,xx); %将xx代⼊解析解，得到解析解对应的数值plot(xx,y,'k--');legend('前向欧拉法','改进欧拉法','欧拉两步法','解析解');oula.m:function f=oula(x,y)f=-3*y+8*x-7;2. 运算结果A1,A2为前向欧拉法在节点处的近似值，B1,B2为改进的欧拉法在节点处的近似值，C1,C2为欧拉公式法在节点处的近似值。
matlab 求解微分方程
在 MATLAB中可以使用 ode45 或者 ode15s 函数来求解常微分方程。
如果想要求解初值问题，可以使用 ode45 函数，语法如下：
```
tspan = [t0, tf]; % t0为初始时刻，tf为结束时刻
y0 = [y1, y2, ..., yn]; % y1, y2, ..., yn为初始条件
[t, y] = ode45(@(t, y) diffeq(t, y), tspan, y0);
```
其中，`diffeq` 是一个用户定义的函数，用来表示微分方程的右端，它的输入参数为时间 t 和状态变量 y，输出为微分方程的右端的值。
`t` 是时间向量，`y` 是状态变量的解。
如果想要求解延迟微分方程或者刚性微分方程，可以使用ode15s 函数，语法和 ode45 函数类似：
```
[t, y] = ode15s(@(t, y) diffeq(t, y), tspan, y0);
```
需要注意的是，求解微分方程之前，需要先定义好微分方程的右端函数 `diffeq` 。
实验七 用matlab 求解常微分方程一、实验目的：1、熟悉常微分方程的求解方法，了解状态方程的概念；2、能熟练使用dsolve 函数求常微分方程（组）的解析解；3、能熟练应用ode45\ode15s 函数分别求常微分方程的非刚性、刚性的数值解；4、掌握绘制相图的方法 二、预备知识： 1．微分方程的概念未知的函数以及它的某些阶的导数连同自变量都由一已知方程联系在一起的方程称为微分方程。
如果未知函数是一元函数，称为常微分方程。
常微分方程的一般形式为0),,",',,()(=n y y y y t F如果未知函数是多元函数，成为偏微分方程。
联系一些未知函数的一组微分方程组称为微分方程组。
微分方程中出现的未知函数的导数的最高阶解数称为微分方程的阶。
若方程中未知函数及其各阶导数都是一次的，称为线性常微分方程，一般表示为)()(')()(1)1(1)(t b y t a y t a y t a y n n n n =++++--若上式中的系数ni t a i ,,2,1),( =均与t 无关，称之为常系数。
2．常微分方程的解析解有些微分方程可直接通过积分求解.例如,一解常系数常微分方程1+=y dt dy可化为dt y dy=+1,两边积分可得通解为1-=tce y .其中c 为任意常数.有些常微分方程可用一些技巧,如分离变量法,积分因子法,常数变异法,降阶法等可化为可积分的方程而求得解析解.线性常微分方程的解满足叠加原理,从而他们的求解可归结为求一个特解和相应齐次微分方程的通解.一阶变系数线性微分方程总可用这一思路求得显式解。
高阶线性常系数微分方程可用特征根法求得相应齐次微分方程的基本解，再用常数变异法求特解。
一阶常微分方程与高阶微分方程可以互化，已给一个n 阶方程),,",',()1()(-=n n y y y t f y设)1(21,,',-===n n y y y y y y ，可将上式化为一阶方程组⎪⎪⎪⎩⎪⎪⎪⎨⎧====-),,,,(''''2113221n n nn y y y t f y yy y y y y反过来，在许多情况下，一阶微分方程组也可化为高阶方程。
用matlab求解常微分方程在MATLAB中，由函数dsolve()解决常微分方程(组)的求解问题，其具体格式如下：r 二dsolve('eql,eq2,・••字condl,cond2,・.・；V)匕ql,eq2,・・・*为微分方程或微分方程组，，condl,cond2,.・・；是初始条件或边界条件，P是独立变量，默认的独立变量是讥函数dsolve用来解符号常微分方程、方程组，如果没有初始条件，则求出通解，如果有初始条件，则求出特解。
dy _1例1：求解常微分方程莎一石的MATLAB程序为：dsolve(* Dy=l/(x+y) 1r!x1),注意，系统缺省的自变量为t,因此这里要把自变量写明。
其中：Y=lambertw(X)表示函数关系Y*exp(Y)二X。
例2：求解常微分方程E'-y— 0的MATLAB程序为：Y2=dsolve(1y*D2y-Dy A2=01, 1x f)Y2=dsolve(!D2y*y-Dy A2=0 J )我们看到有两个解，其中一个是常数0。
dx 心 ？ —+ 5x + y = e dt空_兀_3『= g2f 例3：求常微分方程组〔力 ' 通解的MATLAB 程序为： [X,Y]=dsolve(f Dx+5*x+y=exp(t),Dy-x-3*y=exp(2*t) 1, 111) [X,Y]=dsolve(1Dx+2 *x-Dy=l0 * cos(t),Dx+Dy+2 *y=4 *exp(-2*t) T ,f x(0)=2,y(0)=0f ,f t T) 以上这些都是常微分方程的精确解法，也称为常微分方程的符号解。
但是，我们知 道，有大量的常微分方程虽然从理论上讲，其解是存在的，但我们却无法求出其解析 解，此吋，我们需要寻求方程的数值解，在求常微分方程数值解方面，MATLAB 具有丰 富的函数，我们将其统称为solver,其一般格式为：i°cosr, 7=2 例4：求常微分方程组 y = 0 z 通解的MATLAB 程序为:[T,Y]=solver(odefun,tspan,yO)该函数表示在区间tspan=[tO,tf]±,用初始条件yO求解显式常微分方程卩="，刃。
solver 为命令odc45, odc23, ode 113, ode 15s, ode23s, ode23t, ode23tb 之一,这些命令各有特点。
我们列表说明如下:odefun为显式常微分方程）,”，刃中的“丿）tspan为求解区间，要获得问题在其他指定点心，也，…上的解，则令切期二［心,/“，…心］（要求-单调递增或递减），yO初始条件。
例5：求解常微分方程y' = -2y + 2x2+2x, 0<x<0.5, y（0） = l的MATLAB程序如下：y=dsolve（1Dy=-2*y+2*x A2+2*x1, 1y （0）=11, 1x!） x=0:0.01:0.5;yy=subs(y,x);fun=inline (*-2*y+2*x*x+2*x1); [x,y]=odel5s(fun, [0:0.01:0.5],1);ys=x ・ *x+exp(一2*x);Plot(x,y, f r f,x,ys, f b f )例6：求解常微分方程分冷+尸°"3)"'八°2()的解，并画出解的图形。
分析：这是一个二阶非线性方程(函数以及所有偏导数军委一次幕的是现性方程，高于一次的为非线性方程)，用现成的方法均不能求解，但我们可以通过下面的变换，将二阶方程化为一阶方程组，即可求解。
= dy_令：兀产儿氐-亦，“ =7,则得到：~~ =兀2 ，兀| (0) = 1at<= 7(1-%,2 )x2 - Xj, x2 (0) = 0dt解:function [dfy]=mytt(t,fy)%fl=y;f2=dy/dt%求二阶非线性微分方程时，把一阶、二阶直到(ml)阶导数用另外一个函数代替%用ode45命令时，必须表示成Y』f(t,Y)的形式%Y=[yl;y2;y3],Y,=[yl\y2\y3>[y2;y3;f(yl,y2,y3)],%其中yl=y,y2=y',y3=y"%更高阶时类似dfy=[fy(2);7*(l ・fy ⑴八2)*fy(2)・fy ⑴]；clear;clc[t,yy]=ode45Cmytt\[0 40],[1;0]);plot(t,yy)legend('yTdy')【例4.14.2.1-1］采用ODE解算扌旨令研究围绕地球旋转的卫星轨道。
(1)问题的形成轨道上的卫星，在牛顿第二定律F=ma = m^,和万有引力定律F T叫；作用下有dt r2a^L = -G^r ,引力常数G=6.672*10 H(N.m2/kg2) ,ME=5.97*1024(kg)是地球的质量。
假dr r定卫星以初速度v y(0)=4000m/s在x(0)=4.2*l()7(m)处进入轨道。
(2) 构成一阶微分方程组令 Y 二[刃力力为]「二[尤 V V x v y ]T =[x y * y]T儿 (…严 儿 (宀)计(3) 根据上式[dYdt.m]function Ycl=DYdt (t, Y)% t% Yglobal G ME% xy 二Y(l:2) ;Vxy=Y(3:4) ; %r=sqrt(sum(xy.八2)); Yd 二[Vxy ;-G*ME*xy//3] ; %(4)global G ME % <1>G=6 ・ 672e-ll;ME=5 ・ 97e24;vy0=4000;x0=-4 ・ 2e7;t0=0;tf=60*60*24*9; tspan=[tO,tf];%-GM E —GMY0=[xO；0;0;vyO]; %X=YY(:,1);Y=YY(:,2) ; %plot(X,Y, , linewidth1,2); hold on%axis(* image 1) %[XE,YE,ZE] = sphere (10); %RE=0・64e7; %XE=RE*XE;YE=RE*YE;ZE=0*ZE; %mesh(XE,YE,ZE),hold off %练习：dy_1.利用MATLAB求常微分方程的初值问题杰* " , >?L=o = 2的解。
r=dsolve(1Dy+3*y=01, !y (0)=21, 1x *)2.利用MATLAB求常微分方程的初值问题（1 +扌）* = 23,儿=。
二1,儿乜”的解。
r=dsolve(1D2y*(l+x A2)-2*x*Dy=0',1y(0)=1,Dy(0)=31, *)3. 利用MATLAB 求常微分方程y ⑷-2厂+* = 0的解。
解:y=dsolve(,D4y-2*D3y+D2y ,/x ,) [X,Y]=dsolve(1Dx*2+4*x+Dy-y=exp(t),Dx+3*x+y=01, 1x(0)=1.5,y (0)=01, 丫 t 】) 5.求解常微分方程/，-2(l-/)/+y = 0, 0<x<30, y(0) = l, /(0) = 0的特解，并作出解函 数的曲线图。
r=dsolve(1D2y-2*(l-y A 2)*Dy+y=01, 1y(0)=0,Dy(0)=01, * x 1)4.利用MATLAB 求常微分方程组r dx A dy.2— + 4x + -- y = e dt dt 3 X l /=0 ~ 2 y = 0 /=0的特解。
function DY=mytt2(t,Y)DY=[Y(2);2*( 1 -Y(l )A2)*Y(2)+Y( 1)]; clear;clc[t,YY]=ode45(,mytt2\[0 30],[l;0]); plot(t,YY)legendC^Vdy 1) 求解特殊函数方程勒让德方程的求解dy dx + /(/ + l)y = 0dx d , o x—(1一厂)即dy —2— + /(/ + l)y = 0 dxr=dsolve(1(l-x A2)*D2y-2*x*Dy+y*l*(1+1)=01,1x')连带勒让徳方程的求解:r=dsolve(1(l-x A 2)*D2y-2*x*Dy+y*(1*(l+l)-m A 2/(l-x A 2))=0\)贝塞尔方程r=dsolve(!x A 2*D2y+x*Dy+(x A 2-v A 2)*y=01)利用maplemaple dsolve((l-x A 2)*diff(y(x),x$2)-2*x*diff(y(x)r x)+n*(n+1)*y(x) =0, y (x)) 利用MAPLE 的深层符号计算资源d)等—力孚+ (+1)- 2m i 21-x y = 0经典特殊函数的调用MAPLE库函数在线帮助的检索树发挥MAPLE的计算潜力调用MAPLE函数【例 5.7.3.1-1 ]求递推方程f(n) = -3/(n- 1) -2/(n-2)的通解。
(1)gsl=maple(!rsolve(f(n)=-3*f(n-1)-2*f (n-2),f (k)); 1)(2)调用格式二gs2=maple(T rsolve1, T f(n)=-3*f(n-1)-2*f(n-2) T, !f(k) !)【例5.7.3.1-2]求/* =供的Hessian矩阵。
(1)FHl=maple(!hessian(x*y*z f [x,y,z]); *)FH2=maple(* hessian!, * x*y*z1, 1 [x,y,z] 1)FH=sym(FH2)【例5.7.3」・3】求sin(>2 + y2)在兀二Q y = 0处展开的截断8阶小量的台劳近似式。
(1)TLl=maple(^taylor(sin(x A2+y A2),[x=0,y=0],8)1)(2)maple(1 readlib(mtaylor); 1);TL2=maple(1mtaylor(sin(x A2+y A2), [x=0,y=0],8) 1);pretty(sym(TL2))运行MAPLE程序【例5.7.321】目标：设计求取一般隐函数f(x.y) = 0的导数才⑴解析解的程序，并要求：该程序能象MAPLE原有函数一样可被永久调用。
[DYDZZY.src]DYDZZY:=proc (f)#DYDZZY (f) is used to get the derivate of#an implicit functionlocal Eq,deq,imderiv;Eq: = 1 Eq 1 ;Eq:=f;deq:=diff(Eq, x);readlib(isolate);imderiv:=isolate(deq,diff(y(x) , x));end;procread(1DYDZZY•src1)(3)sl=maple(1DYDZZY(x=log(x+y(x)));1)s2=maple(f DYDZZY(x A2*y(x)-exp(2*x)=sin(y(x))) !) s3=maple(1DYDZZY1, 1 cos(x+sin(y(x)))=sin(y(x)) 1)clear maplemexprocread(1DYDZZY ・ src1);maple(1 save(、DYDZZY.m s) 1);(5)maple(* read'J'DYDZZY ・m、T);ss2=maple(1DYDZZY(x A2*y(x)-exp(2*x)=sin(y(x)))1)。
