(完整版)实验七用matlab求解常微分方程
matlab-常微分方程

Events
含义 为‘on’时,控制解向量 有效值: 范数的相对误差,使每 on、off 步计算中,满足: 缺省值: norm(e)<=max(RelTol*n off orm(y),AbsTol) 有效值: 为‘on’时,返回相应的 on、off 事件记录
取值
参数设置
属性名 含义 若无输出参量,则solver 将执行下面操作之一: 有效值: 画出解向量中各元素随 odeplot、 时间的变化; odephas2、画出解向量中前两个分 odephas3、量构成的相平面图; odeprint 画出解向量中前三个分 缺省值: 量构成的三维相空间图 odeplot ; 随计算过程,显示解向 量 取值
使用于精度较低 的情形
OD 求解器 E类 Solver 型 非 ode113 刚 性 适 度 ode23t 刚 性 刚 ode15s 性
特点
说明
多步法;Adams算 计算时间比ode45 法;高低精度均可 短 -3~10-6 到10 采用梯形算法 适度刚性情形
多步法;Gear’s反 若ode45失效时, 向数值微分;精度 可尝试使用 中等
• (6)若没有给定输出参量,则在命令窗口显 示解列表。若该命令找不到解析解,则返 回一警告信息,同时返回一空的sym对象。 这时,用户可以用命令ode23或ode45求解 方程组的数值解。
y ′′ = −a y ′ y ( 0) = 1 y ′(π / a ) = 0
2
例1
例2
• >> [u,v] = dsolve('Du=v,Dv=u') u= C1*exp(-t)+C2*exp(t) V= -C1*exp(-t)+C2*exp(t)
常微分方程(ODEs)的MATLAB数值解法

1.ODE 解算器简介 ..............................................................................................................................................................3 2.微分方程转换 ...................................................................................................................................................................5 3.刚性/非刚性问题 ..............................................................................................................................................................8 4.隐式微分方程(IDE) ........................................................................................................................................................10 5.微分代数方程(DAE)....................................................................................................................................................... 15 6.延迟微分方程(DDE)....................................................................................................................................................... 18 7.边值问题(BVP) ...............................................................................................................................................................20
matlab常微分方程的数值解法实验报告

实验四常微分方程的数值解法 指令:[t,y]=ode23(‘fun ’,tspan,yo) 2/3阶龙格库塔方法 [t,y]=ode45(‘fun ’,tspan,yo) 4/5阶龙格库塔方法 [t,y]=ode113(‘fun ’,tspan,yo) 高阶微分方程数值方法其中fun 是定义函数的文件名。
该函数fun 必须以为dx 输出量,以t,y 为输入量。
tspan=[t0 tfina]表示积分的起始值和终止值。
yo 是初始状态列向量。
考虑到初始条件有00d , (0)0,d d , (0)0.d SSI S S tI SI I I I tββμ⎧=-=>⎪⎪⎨⎪=-=≥⎪⎩ (5.24) 这就是Kermack 与McKendrick 的SIR 仓室模型. 方程(5.24)无法求出()S t 和()I t 的解析解.我们先做数值计算。
Matlab 代码为:function dy=rigid(t,y) dy=zeros(2,1); a=1; b=0.3;dy(1)=a*y(1).*y(2)-b*y(1); dy(2)=-a*y(1).*y(2);ts=0:.5:50; x0=[0.02,0.98];[T,Y]=ode45('rigid',ts,x0); %plot(T,Y(:,1),'-',T,Y(:,2),'*') plot(Y(:,2),Y(:,1),'b--') xlabel('s') ylabel('i')任务:1 画出i (t ),2分析各参数的影响例57:求解两点边值问题:0)5(,0)1(,32==='-''y y x y y x 。
(注意:相应的数值解法比较复杂)。
y=dsolve('x*D2y-3*Dy=x^2','y(1)=0,y(5)=0','x') ↙ y =-1/3*x^3+125/468+31/468*x^4例:用数值积分的方法求解下列微分方程 π21''2t y y -=+设初始时间t0=0;终止时间tf=3*pi ;初始条件0|',0|00====x x y y 。
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程序以下是一个简单的MATLAB 程序,用于求解一阶常微分方程:matlab复制代码% 定义微分方程 dy/dx = f(x, y)f = @(x, y) -x*y;% 初始条件 y(0) = 1y0 = 1;% 定义 x 的范围xspan = [0, 10];% 使用 MATLAB 内置函数 ode45 进行求解[t, y] = ode45(f, xspan, y0);% 绘制解的图形plot(t, y(:,1));xlabel('x');ylabel('y');title('Solution of the differential equation dy/dx = -xy');在这个程序中,我们定义了一个一阶常微分方程dy/dx = -xy,并使用MATLAB 内置函数ode45进行求解。
初始条件为y(0) = 1,求解范围为xspan = [0, 10]。
最后,我们使用plot函数绘制了解的图形。
这个程序是用来求解一阶常微分方程的,而这个方程是dy/dx = -xy。
这是一个简单的线性方程,但它的解在物理和工程中有许多实际应用。
接下来,我们逐行解释一下代码:1.% 定义微分方程 dy/dx = f(x, y):这是一个注释,说明下面的代码是定义微分方程。
2. f = @(x, y) -x*y;:这行定义了一个匿名函数f,它接受两个参数x和y,并返回-x*y。
这个函数就是我们的微分方程dy/dx的右边部分。
3.% 初始条件 y(0) = 1:这是一个注释,说明下面的代码是定义初始条件。
4.y0 = 1;:这行定义了初始条件y(0) = 1,也就是说当x=0时,y=1。
5.% 定义 x 的范围:这是一个注释,说明下面的代码是定义自变量x的范围。
6.xspan = [0, 10];:这行定义了自变量x的范围从0到10。
7.% 使用 MATLAB 内置函数 ode45 进行求解:这是一个注释,说明下面的代码将使用MATLAB 的内置函数ode45来求解微分方程。
matlab求解常微分方程

matlab求解常微分⽅程本⽂主要介绍matlab中求解常微分⽅程(组)的dsolve和ode系列函数,并通过例⼦加深读者的理解。
⼀、符号介绍D: 微分符号;D2表⽰⼆阶微分,D3表⽰三阶微分,以此类推。
⼆、函数功能介绍及例程1、dsolve 函数dsolve函数⽤于求常微分⽅程组的精确解,也称为常微分⽅程的符号解。
如果没有初始条件或边界条件,则求出通解;如果有,则求出特解。
1)函数格式Y = dsolve(‘eq1,eq2,…’ , ’cond1,cond2,…’ , ’Name’)其中,‘eq1,eq2,…’:表⽰微分⽅程或微分⽅程组;’cond1,cond2,…’:表⽰初始条件或边界条件;‘Name’:表⽰变量。
没有指定变量时,matlab默认的变量为t;2)例程例1.1(dsolve 求解微分⽅程)求解微分⽅程:dsolve('Dy=3*x^2','x')例1.2(加上初始条件)求解微分⽅程:例2(dsolve 求解微分⽅程组)求解微分⽅程组:由于x,y均为t的导数,所以不需要在末尾添加’t’。
2、ode函数在上⽂中我们介绍了dsolve函数。
但有⼤量的常微分⽅程,虽然从理论上讲,其解是存在的,但我们却⽆法求出其解析解,此时,我们需要寻求⽅程的数值解。
ode是Matlab专门⽤于解微分⽅程的功能函数。
该求解器有变步长(variable-step)和定步长(fixed-step)两种类型。
不同类型有着不同的求解器,具体说明如下图。
其中,ode45求解器属于变步长的⼀种,采⽤Runge-Kutta算法;其他采⽤相同算法的变步长求解器还有ode23。
ode45表⽰采⽤四阶-五阶Runge-Kutta算法,它⽤4阶⽅法提供候选解,5阶⽅法控制误差,是⼀种⾃适应步长(变步长)的常微分⽅程数值解法,其整体截断误差为(Δx)^5。
解决的是Nonstiff(⾮刚性)常微分⽅程。
matlab求解常微分方程.docx

用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求解显式常微分方程卩=",刃。
实验报告七常微分方程初值问题的数值解法

浙江大学城市学院实验报告课程名称数值计算方法实验项目名称常微分方程初值问题的数值解法 实验成绩指导老师签名日期2015/12/16 一.实验目的和要求1. 用Matlab 软件掌握求微分方程数值解的欧拉方法和龙格-库塔方法; 2. 通过实例学习用微分方程模型解决简化的实际问题;二.实验内容和原理编程题2-1要求写出Matlab 源程序m 文件,并有适当的注释语句;分析应用题2-2,2-3,2-4,2-5要求将问题的分析过程、Matlab 源程序和运行结果和结果的解释、算法的分析写在实验报告上; 2-1 编程编写用向前欧拉公式和改进欧拉公式求微分方程数值解的Matlab 程序,问题如下:在区间[],a b 内(1)N +个等距点处,逼近下列初值问题的解,并对程序的每一句添上注释语句; Euler 法y=eulera,b,n,y0,f,f1,b1改进Euler 法y=eulerproa,b,n,y0,f,f1,b1 2-2 分析应用题假设等分区间数100n =,用欧拉法和改进欧拉法在区间[0,10]t ∈内求解初值问题()()20(0)10y t y t y '=-⎧⎨=⎩并作出解的曲线图形,同时将方程的解析解也画在同一张图上,并作比较,分析这两种方法的精度; 2-3 分析应用题用以下三种不同的方法求下述微分方程的数值解,取10h = 画出解的图形,与精确值比较并进行分析; 1欧拉法; 2改进欧拉法; 3龙格-库塔方法;2-4 分析应用题考虑一个涉及到社会上与众不同的人的繁衍问题模型;假设在时刻t 单位为年,社会上有人口()x t 人,又假设所有与众不同的人与别的与众不同的人结婚后所生后代也是与众不同的人;而固定比例为r 的所有其他的后代也是与众不同的人;如果对所有人来说出生率假定为常数b ,又如果普通的人和与众不同的人的婚配是任意的,则此问题可以用微分方程表示为:其中变量()()()i p t x t x t =表示在时刻t 社会上与众不同的人的比例,()i x t 表示在时刻t 人口中与众不同的人的数量;1假定(0)0.01,0.02p b ==和0.1r =,当步长为1h =年时,求从0t =到50t =解()p t 的近似值,并作出近似解的曲线图形;2精确求出微分方程的解()p t ,并将你当50t =时在分题b 中得到的结果与此时的精确值进行比较; MATLAB 相关函数求微分方程的解析解及其数值的代入dsolve‘egn1’,‘egn2’,‘x ’ subsexpr,{x,y,…},{x1,y1,…}其中‘egn i ’表示第i 个方程,‘x ’表示微分方程中的自变量,默认时自变量为t ; subs 命令中的expr 、x 、y 为符合型表达式,x 、y 分别用数值x1、x2代入; >>symsxyz>>subs'x+y+z',{x,y,z},{1,2,3} ans= 6>>symsx>>subs'x^2',x,2 ans= 4>>s=dsolve‘12Dy y ∧=+’,‘(0)1y =’,‘x ’ ans= >>symsx >>subss,x,2 ans=右端函数(,)f x y 的自动生成f=inline ‘expr ’,’var1’,‘var2’,……其中’expr ’表示函数的表达式,’var1’,‘var2’表示函数表达式中的变量,运行该函数,生成一个新的函数表达式为fvar1,var2,……; >>f=inline'x+3y','x','y' f=Inlinefunction: fx,y=x+3y >>f2,3 ans= 114,5阶龙格-库塔方法求解微分方程数值解t,x=ode45f,ts,x0,options其中f 是由待解方程写成的m 文件名;x0为函数的初值;t,x 分别为输出的自变量和函数值列向量,t的步长是程序根据误差限自动选定的;若ts=t0,t1,t2,…,tf,则输出在自变量指定值,等步长时用ts=t0:k:tf,输出在等分点;options 用于设定误差限可以缺省,缺省时设定为相对误差310-,绝对误差610-,程序为:options=odeset ‘reltol ’,rt,’abstol ’,at,这里rt,at 分别为设定的相对误差和绝对误差;常用选项见下表;选项名 功能 可选值 省缺值 AbsTol 设定绝对误差正数 RelTol 设定相对误差 正数InitialStep 设定初始步长 正数 自动 MaxStep设定步长上界正数MaxOrder 设定ode15s 的最高阶数 1,2,3,4,5 5 Stats 显示计算成本统计 on,off off BDF 设定ode15s 是否用反向差分on,offoff例:在命令窗口执行>>odefun =inline ‘2*y t y -’,‘t ’,‘y ’;>>[],45(,[0,4],1)t y ode odefun =;ans=>>t y ‘o-’,%解函数图形表示>>45(,[0,4],1)ode odefun %不用输出变量,则直接输出图形 >>[],45(,0:4,1)t y ode odefun =;[],t yans=三.操作方法与实验步骤包括实验数据记录和处理2-1编程编写用向前欧拉公式和改进欧拉公式求微分方程数值解的Matlab 程序,问题如下:在区间[],a b 内(1)N +个等距点处,逼近下列初值问题的解,并对程序的每一句添上注释语句; Euler 法y=eulera,b,n,y0,f,f1,b1改进Euler 法y=eulerproa,b,n,y0,f,f1,b1Euler 法y=eulera,b,n,y0,f,f1,b1 y=zeros1,n+1; y1=y0; h=b-a/n; x=a:h:b; fori=1:n; yi+1=yi+hfxi,yi; end plotx,y holdon%求微分方程的精确解 x1=linspacea,b,100; '精确解为' s=dsolvef1,b1,'x' symsxy1=zeros1,100; for i=1:100y1i=subss,x,x1i; endplotx1,y1,'r'title'红色代表精确解'改进Euler 法y=eulerproa,b,n,y0,f,f1,b1 %求微分方程的数值解 y=zeros1,n+1; y1=y0; h=b-a/n; x=a:h:b; fori=1:n; T1=fxi,yi; T2=fxi+1,yi+hT1; yi+1=yi+h/2T1+T2; end plotx,y holdon%求微分方程的精确解 x1=linspacea,b,100; '精确解为' s=dsolvef1,b1,'x' symsxy1=zeros1,100; fori=1:100 y1i=subss,x,x1i; endplotx1,y1,'r'title'红色代表精确解' 2-2分析应用题假设等分区间数100n =,用欧拉法和改进欧拉法在区间[0,10]t ∈内求解初值问题()()20(0)10y t y t y '=-⎧⎨=⎩并作出解的曲线图形,同时将方程的解析解也画在同一张图上,并作比较,分析这两种方法的精度;1向前欧拉法>>euler0,10,100,10,inline'y-20','x','y','Dy=y-20','y0=10' ans= 精确解为 s= 20-10expx ans= +005Columns1through8(2)改进欧拉法>>eulerpro0,10,100,10,inline'y-20','x','y','Dy=y-20','y0=10' ans= 精确解为 s= 20-10expx ans= +005Columns1through8改进欧拉法的精度比向前欧拉法更高; 2-3分析应用题用以下三种不同的方法求下述微分方程的数值解,取10h = 画出解的图形,与精确值比较并进行分析; 1欧拉法; 2改进欧拉法;2-4分析应用题考虑一个涉及到社会上与众不同的人的繁衍问题模型;假设在时刻t 单位为年,社会上有人口()x t 人,又假设所有与众不同的人与别的与众不同的人结婚后所生后代也是与众不同的人;而固定比例为r 的所有其他的后代也是与众不同的人;如果对所有人来说出生率假定为常数b ,又如果普通的人和与众不同的人的婚配是任意的,则此问题可以用微分方程表示为:其中变量()()()i p t x t x t =表示在时刻t 社会上与众不同的人的比例,()i x t 表示在时刻t 人口中与众不同的人的数量;1假定(0)0.01,0.02p b ==和0.1r =,当步长为1h =年时,求从0t =到50t =解()p t 的近似值,并作出近似解的曲线图形;2精确求出微分方程的解()p t ,并将你当50t =时在分题b 中得到的结果与此时的精确值进行比较;1>>euler0,50,50,,inline'','t','p','Dp=','p0= 1' ans= 精确解为 s=1-99/100expx/500 ans=Columns1through82>>dsolve'Dp=','p0=','t' ans=1-99/100expt/500 >>1-99/100exp ans=与欧拉法求得的精确值差0,0001四.实验结果与分析。
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
实验七 用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反过来,在许多情况下,一阶微分方程组也可化为高阶方程。
所以一阶微分方程组与高阶常微分方程的理论与方法在许多方面是相通的,一阶常系数线性微分方程组也可用特征根法求解。
3.微分方程的数值解法除常系数线性微分方程可用特征根法求解,少数特殊方程可用初等积分法求解外,大部分微分方程无限世界,应用中主要依靠数值解法。
考虑一阶常微分方程初值问题⎩⎨⎧=<<=000)()),(,()('y t y t t t t y t f t y f其中)'.,,,(,)',,,(,)',,,(020*******m m m y y y y f f f f y y y y ===所谓数值解法,就是寻求)(t y 在一系列离散节点f n t t t t ≤<<< 10上的近似值nk y k ,,1,0, =称kk k t t h -=+1为步长,通常取为常量h 。
最简单的数值解法是Euler 法。
Euler 法的思路极其简单:在节点出用差商近似代替导数h t y t y t y k k k )()()('1-≈+这样导出计算公式(称为Euler 格式),2,1,0),,(1=+=+k y t hf y y k k k k他能求解各种形式的微分方程。
Euler 法也称折线法。
Euler 方法只有一阶精度,改进方法有二阶Runge-Kutta 法、四阶Runge-Kutta 法、五阶Runge-Kutta-Felhberg 法和先行多步法等,这些方法可用于解高阶常微分方程(组)初值问题。
边值问题采用不同方法,如差分法、有限元法等。
数值算法的主要缺点是它缺乏物理理解。
4.解微分方程的MATLAB 命令MATLAB 中主要用dsolve 求符号解析解,ode45,ode23,ode15s 求数值解。
ode45类似,只是精度低一些。
ode12s 用来求解刚性方程组,是用格式同ode45。
可以用help dsolve, help ode45查阅有关这些命令的详细信息. 例1 求下列微分方程的解析解(1)b ay y +='(2)1)0(',0)0(,)2sin(''==-=y y y x y(3)1)0(',1)0(',','==-=+=gffgggff方程(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) %以最简形式显示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)例2求解微分方程,1)0(,1'=++-=ytyy先求解析解,再求数值解,并进行比较。
由>>clear;>>s=dsolve('Dy=-y+t+1','y(0)=1','t')>>simplify(s)可得解析解为tety-+=。
下面再求其数值解,先编写M文件fun8.m%M函数fun8.mfunction 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')结果见图7.1图16.6.1 解析解与数值解由图16.6.1可见,解析解和数值解吻合得很好。
例3 求方程)0(',)0(,sin "0===θθθθθmg ml的数值解.不妨取15)0(,8.9,1===θg l .则上面方程可化为0)0(',15)0(,sin 8.9"===θθθθ先看看有没有解析解.运行MATLAB 代码>>clear;>>s=dsolve('D2y=9.8*sin(y)','y(0)=15','Dy(0)=0','t') >>simplify(s)知原方程没有解析解.下面求数值解.令',21θθ==y y 可将原方程化为如下方程组⎪⎩⎪⎨⎧====0)0(,15)0()sin(8.9''211221y y y y y y建立M 文件fun9.m 如下%M 文件fun9.mfunction 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')结果见图7.2图7.2 数值解由图7.2可见,θ随时间t周期变化。
以上这些都是常微分方程的精确解法,也称为常微分方程的符号解。
但是,我们知道,有大量的常微分方程虽然从理论上讲,其解是存在的,但我们却无法求出其解析解,此时,我们需要寻求方程的数值解,在求常微分方程数值解方面,MATLAB具有丰富的函数,我们将其统称为solver,其一般格式为:[T,Y]=solver(odefun,tspan,y0)该函数表示在区间tspan=[t0,tf]上,用初始条件y0求解显式常微分方程y f t y='(,)solver为命令ode45,ode23,ode113,ode15s,ode23s,ode23t,ode23tb之一,这些命令各有特点。
我们列表说明如下:odefun 为显式常微分方程'(,)y f t y =中的(,)f t y tspan 为求解区间,要获得问题在其他指定点012,,,t t t 上的解,则令012[,,,,]f tspan t t t t =(要求i t 单调),y0初始条件。
例5:求解常微分方程2'222y y x x =-++,00.5x ≤≤,(0)1y =的MATLAB 程序如下:fun=inline('-2*y+2*x*x+2*x');[x,y]=ode23(fun,[0,0.5],1) 结果为: x =0,0.0400,0.0900,0.1400,0.1900,0.2400,0.2900,0.3400,0.3900,0.4400,0.4900,0.5000 y =1.0000,0.9247,0.8434,0.7754,0.7199,0.6764,0.6440,0.6222,0.6105,0.6084,0.6154,0.6179例6:求解常微分方程222(1)0,(0)1,'(0)0d y dy y y y y dt dt μ--+===的解,并画出解的图形。