(整理)实验五用matlab求解常微分方程.
matlab求解矩阵常微分方程
matlab求解矩阵常微分方程标题：基于Matlab的矩阵常微分方程求解导言：矩阵常微分方程是指微分方程中的未知函数为矩阵，常见于控制系统、信号处理、计算机图形学等领域。
Matlab作为一款功能强大的数学软件，提供了许多工具和函数用于求解矩阵常微分方程。
本文将介绍如何使用Matlab求解矩阵常微分方程的基本方法和步骤。
一、矩阵常微分方程的定义矩阵常微分方程可以写成如下形式：$$\frac{dX(t)}{dt} = A(t)X(t)$$其中，$X(t)$是未知函数矩阵，$A(t)$是已知函数矩阵。
二、Matlab中的矩阵常微分方程求解函数Matlab提供了几个函数用于求解矩阵常微分方程，其中最常用的是ode45函数。
ode45函数是基于龙格-库塔法的数值求解器，可以高效地求解大多数常微分方程问题。
三、使用ode45函数求解矩阵常微分方程的步骤1. 定义矩阵常微分方程的参数和初始条件。
2. 编写一个函数，用于描述矩阵常微分方程的右端项。
3. 调用ode45函数，传入函数句柄和初始条件，求解矩阵常微分方程。
四、示例：求解矩阵常微分方程假设我们要求解如下矩阵常微分方程：$$\frac{dX(t)}{dt} = \begin{bmatrix} 0 & 1 \\ -1 & 0 \end{bmatrix} X(t)$$初始条件为$X(0) = \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix}$。
根据上述步骤，我们可以先定义矩阵常微分方程的参数和初始条件：```matlabA = [0 1; -1 0];X0 = [1 0; 0 1];```然后编写一个函数，用于描述矩阵常微分方程的右端项：```matlabfunction dXdt = matrixODE(t, X)A = [0 1; -1 0];dXdt = A * X;end```最后调用ode45函数求解矩阵常微分方程：```matlab[t, X] = ode45(@matrixODE, [0 10], X0);```其中，@matrixODE表示函数句柄，[0 10]表示求解时间范围。
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 程序，用于求解一阶常微分方程：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梯形法？它又是如何应用于求解常微分方程的呢？让我们来深入了解matlab梯形法。
在matlab中，梯形法是一种常用的数值求解方法，它可以用于求解常微分方程。
该方法的基本思想是将微分方程中的导数用差分代替，从而将微分方程转化为代数方程组，再利用matlab进行求解。
通过该方法，我们可以得到微分方程的数值解，从而更好地理解和分析问题。
现在，让我们来探讨matlab梯形法在求解常微分方程中的应用。
假设我们需要求解如下的一阶常微分方程：\[ \frac{dy}{dt} = f(t,y) \]其中，\( f(t,y) \) 是关于\( t \)和\( y \)的函数。
我们需要将微分方程离散化，即用差分代替导数。
通过将时间区间\( [a, b] \)进行均匀划分，我们可以得到：\[ t_0 = a, t_1, t_2, ..., t_n = b \]\[ y_0 = \alpha, y_1, y_2, ..., y_n \]\[ h = \frac{b-a}{n} \]其中，\( t_i \) 是时间节点，\( y_i \) 是对应的近似解，\( h \) 是时间步长。
接下来，我们可以利用梯形法进行求解。
梯形法的迭代公式为：\[ y_{i+1} = y_i + \frac{h}{2}[f(t_i, y_i) + f(t_{i+1}, y_{i+1})] \]通过不断迭代，我们可以得到微分方程的数值解。
在实际应用中，matlab梯形法可以很好地处理各种类型的常微分方程。
无论是线性方程还是非线性方程，matlab梯形法都能提供较为准确的数值解。
该方法还可以用于求解初值问题和边值问题，具有较好的通用性和适用性。
MATLAB求解常微分方程数值解
利用MATLAB求解常微分方程数值解目录1.容简介12.Euler Method(欧拉法)求解12.1.显式Euler法和隐式Euler法22.2.梯形公式和改进Euler法32.3.Euler法实用性53.Runge-Kutta Method(龙格库塔法)求解53.1.Runge-Kutta基本原理63.2.MATLAB中使用Runge-Kutta法的函数74.使用MATLAB求解常微分方程84.1.使用ode45函数求解非刚性常微分方程84.2.刚性常微分方程95.总结9参考文献11附录121.显式Euler法数值求解122.改进Euler法数值求解123.四阶四级Runge-Kutta法数值求解134.使用ode45求解141.容简介把《高等工程数学》看了一遍，增加对数学容的了解，对其中数值解法比较感兴趣，这大概是因为在其它各方面的学习和研究中经常会遇到数值解法的问题。
理解模型然后列出微分方程，却对着方程无从下手，无法得出精确结果实在是让人难受的一件事情。
实际问题中更多遇到的是利用数值法求解偏微分方程问题，但考虑到先从常微分方程下手更为简单有效率，所以本文只研究常微分方程的数值解法。
把一个工程实际问题弄出精确结果远比弄清楚各种细枝末节更有意思，因此文章中不追求非常严格地证明，而是偏向如何利用工具实际求解出常微分方程的数值解，力求将课程上所学的知识真正地运用到实际方程的求解中去，在以后遇到微分方程的时候能够熟练运用MATLAB得到能够在工程上运用的结果。
文中求解过程中用到MATLAB进行数值求解，主要目的是弄清楚各个函数本质上是如何对常微分方程进行求解的，对各种方法进行MATLAB编程求解，并将求得的数值解与精确解对比，其中源程序在附录中。
最后考察MATLAB中各个函数的适用围，当遇到实际工程问题时能够正确地得到问题的数值解。
2.Euler Method(欧拉法)求解Euler法求解常微分方程主要包括3种形式，即显式Euler法、隐式Euler法、梯形公式法，本节容分别介绍这3种方法的具体容，并在最后对3种方法精度进行对比，讨论Euler法的实用性。
matlab数值求解常微分方程快速方法
MATLAB是一种用于科学计算和工程应用的高级编程语言和交互式环境。
它在数学建模、模拟和分析等方面有着广泛的应用。
在MATLAB 中，常微分方程的数值求解是一个常见的应用场景。
在实际工程问题中，通常需要对常微分方程进行数值求解来模拟系统的动态行为。
本文将介绍MATLAB中对常微分方程进行数值求解的快速方法。
1. 基本概念在MATLAB中，可以使用ode45函数来对常微分方程进行数值求解。
ode45是一种常用的Runge-Kutta法，它可以自适应地选取步长，并且具有较高的数值精度。
使用ode45函数可以方便地对各种类型的常微分方程进行求解，包括一阶、高阶、常系数和变系数的微分方程。
2. 函数调用要使用ode45函数进行常微分方程的数值求解，需要按照以下格式进行函数调用：[t, y] = ode45(odefun, tspan, y0)其中，odefun表示用于描述微分方程的函数，tspan表示求解的时间跨度，y0表示初值条件，t和y分别表示求解得到的时间序列和对应的解向量。
3. 示例演示为了更好地理解如何使用ode45函数进行常微分方程的数值求解，下面我们以一个具体的例子来进行演示。
考虑如下的一阶常微分方程：dy/dt = -2*y其中，y(0) = 1。
我们可以编写一个描述微分方程的函数odefun：function dydt = odefun(t, y)dydt = -2*y;按照上述的函数调用格式，使用ode45函数进行求解：tspan = [0 10];y0 = 1;[t, y] = ode45(odefun, tspan, y0);绘制出解曲线：plot(t, y);4. 高级用法除了基本的函数调用方式外，MATLAB中还提供了更多高级的方法来对常微分方程进行数值求解。
可以通过设定选项参数来控制数值求解的精度和稳定性，并且还可以对刚性微分方程进行求解。
5. 性能优化在实际工程应用中，常常需要对大规模的常微分方程进行数值求解。
Matlab求解常微分方程
Matlab求解常微分方程如果未知函数是一元的，那么对应的微分方程称为常微分方程。
常微分方程的求解问题分为：初值问题和边值问题，以及求解微分代数方程组问题。
一般情况下，只有很少一类常微分方程可以得到解析解。
这里只介绍常微分方程的数值解问题，即根据常微分方程的不同性质，求其离散化的解。
部分常微分方程命令ode23：解非刚性微分方程，低精度，使用Runge-Kutta（龙格-库塔）法的二三阶算法。
ode45：解非刚性微分方程，中等精度，使用Runge-Kutta（龙格-库塔）法的四五阶算法。
ode113：解非刚性微分方程，变精度变阶次Adams-Bashforth-MoultonPECE（阿姆斯特）算法。
ode23t：解中等刚性微分方程，使用自由内插法的梯形法则（多步长外推法）。
ode15s：解刚性微分方程，使用可变阶次的数值微分（NDFs）算法。
ode23s：解刚性微分方程，低阶方法，使用修正的Rosenbrock（罗森伯里克）公式。
ode23tb：解刚性微分方程，低阶方法，使用TR-BDF2方法（方向内推法），即Runger-Kutta公式的第一级采用梯形法则，第二级采用Gear法。
解常微分方程命令的一般语法：[t,y]=ode*('odefun',tspan,y0)或[t,y]=ode*(@odefun,tspan,yo)[t,y]=ode*('odefun',tspan,y0,options)或[t,y]=ode*(@odefun,tspan,yo,options)[t,y]=ode*('odefun',tspan,y0,options,p1,p2,...)[t,y,Te,Ye,Ie]=ode*('odefun',tspan,yo,options)这里ode*：常微分方程命令odefun：常微分方程定义函数tspan：时域t的范围y0：初始值options：解微分方程的各种选项p1,p2,...：微分方程odefun的参数赋值刚性方程：刚性是指其Jacobian矩阵的特征值相差十分悬殊。
用MATLAB求解微分方程
1. 微分方程的解析解
求微分方程（组）的解析解命令:
dsolve(‘方程1’, ‘方程2’,…‘方程n’, ‘初始条件’, ‘自变量’)
结 果：u = tan(t-c)
解 输入命令：dsolve('Du=1+u^2','t')
STEP2
STEP1
解 输入命令: y=dsolve('D2y+4*Dy+29*y=0','y(0)=0,Dy(0)=15','x')
导弹追踪问题
设位于坐标原点的甲舰向位于x轴上点A(1, 0)处的乙舰发射导弹，导弹头始终对准乙舰.如果乙舰以最大的速度v0(是常数)沿平行于y轴的直线行驶，导弹的速度是5v0，求导弹运行的曲线方程.又乙舰行驶多远时，导弹将它击中？ 解法一（解析法）
由(1),(2)消去t整理得模型:
解法二(数值解)
结 果 为：x = (c1-c2+c3+c2e -3t-c3e-3t)e2t y = -c1e-4t+c2e-4t+c2e-3t-c3e-3t+c1-c2+c3)e2t z = (-c1e-4t+c2e-4t+c1-c2+c3)e2t
2、取t0=0，tf=12，输入命令： [T,Y]=ode45('rigid',[0 12],[0 1 1]); plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+')
3、结果如图
图中，y1的图形为实线，y2的图形为“*”线，y3的图形为“+”线.
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
实验五 用matlab 求解常微分方程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是最常用的求解微分方程数值解的命令，对于刚性方程组不宜采用。
ode23与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(',','==-=+=g f f g g g f f 方程（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)例２求解微分方程,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 周期变化。
习题1．求下列微分方程的解析解2．求方程3)0(',1)0(,'2")1(2===+y y xy y x的解析解和数值解,并进行比较3．分别用ode45和ode15s 求解Van-del-Pol 方程()⎪⎩⎪⎨⎧===---1)0',0)0(0)1(1000222x x x dt dxx dtx d的数值解,并进行比较.。