常微分方程(ODEs)的MATLAB数值解法
重要：MATLAB常微分方程(组)数值解法
Matlab常微分方程求解问题分类
边值问题：
初值问题：
• 定解附加条件在自变量 的一端
• 一般形式为： y' f (x, y)
y(a)
y0
• 初值问题的数值解法一 般采用步进法，如 Runge-Kutta法
➢ 在自变量两端均给定附加 条件
y' f (x, y)
➢ 一般形式：y(a)y1, y(b)y2
1.根据常微分方程要求的求解精度与速度要求
求解初值问题：
y
'
y
2x y
y ( 0 ) 1
(0x1)
比较ode45和ode23的求解精度和速度
ode45和ode23的比较-1
function xODE clear all clc
format long
y0 = 1; [x1,y1] = ode45(@f,[0,1],y0); [x2,y2] = ode23(@f,[0,1],y0); plot(x1,y1,'k-',x2,y2,'b--') xlabel('x') ylabel('y')
rD = k(3)*C(2)-k(5)*C(4);
rE = k(4)*C(3)+k(5)*C(4);
% Mass balances dCdt = [rA; rB; rC; rD; rE];
三个串联的CSTR等温反应器（例4-3）
function IsothermCSTRs clear all clc CA0 = 1.8; % kmol/m^3 CA10 = 0.4; % kmol/m^3 CA20 = 0.2; % kmol/m^3 CA30 = 0.1; % kmol/m^3 k = 0.5; % 1/min tau = 2; stoptime = 2.9; % min [t,y] = ode45(@Equations,[0 stoptime],[CA10 CA20 CA30],[],k,CA0,tau); disp(' Results:') disp(' t CA1 CA2 CA3') disp([t,y]) plot(t,y(:,1),'k--',t,y(:,2),'b:',t,y(:,3),'r-') legend('CA_1','CA_2','CA_3') xlabel('Time (min)') ylabel('Concentration') % -----------------------------------------------------------------function dydt = Equations(t,y,k,CA0,tau) CA1 = y(1); CA2 = y(2); CA3 = y(3); dCA1dt = (CA0-CA1)/tau - k*CA1; dCA2dt = (CA1-CA2)/tau - k*CA2; dCA3dt = (CA2-CA3)/tau - k*CA3; dydt = [dCA1dt; dCA2dt; dCA3dt];
常微分方程(ODEs)的MATLAB数值解法
1.ODE 解算器简介 ..............................................................................................................................................................3 2.微分方程转换 ...................................................................................................................................................................5 3.刚性/非刚性问题 ..............................................................................................................................................................8 4.隐式微分方程(IDE) ........................................................................................................................................................10 5.微分代数方程(DAE)....................................................................................................................................................... 15 6.延迟微分方程(DDE)....................................................................................................................................................... 18 7.边值问题(BVP) ...............................................................................................................................................................20
matlab常微分方程数值解法精品PPT课件
2012
第七讲 常微分方程数值解法
内容提要
引言 欧拉近似方法 龙格-库塔(R-K)方法 MATLAB的常微分方程函数 小结
1、引言
物理学所研究的各种物质运动中，有许多物质运动的 过程是用常微分方程来描述的。
例如,质点的加速运动，简谐振动等。
F m dv dt
d2x 2x2 0
f xn1
x0
y(x), x dx
y(xn ) f y(xn ), xn h
作如下近似
yn y(xn )
得：
yn1 yn f yn , xn h
2.1.4 欧拉法误差
利用泰勒级数得：
y xn1 y(xn h)
y(xn )
hy(xn )
1 2
h2
y(xn )
y(xn )
x2 …. xn ….
y(x0) y(x1) y(x2) …. y(xn) ….
y0
y1 y2 …. yn ….
在xn节点上，微分方程可以写为
y(xn1) y(xn ) f y(xn ) , xn h
作如下近似：
yn y(tn )
则得到欧拉解法递推公式的一般形式：
yn1 yn f ( yn , xn ) h
hf
y(xn ),
xn
1 2
h2 y(xn )
作如下近似
yn y(xn )
yn1 yn f yn , xn h
局部截断误差
y
y0
dy dx
x0 x x0
y0 f ( y0 , x0 ) x x0
此切线与x=x1交点纵坐标为：
y1 y0 f ( y0 , x0 ) x1 x0
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_常微分方程数值解法
dt 2
简朴问题可以求得解析解，多数实际问题靠数值求解 。
第4页
一阶常微分方程(ODE )初值问题 ： ODE ：Ordinary Differential Equation
dy
f
(x,
y)
dx
x0 x xn
y(x0 ) y0
数值解法就是求y(x)在某些分立旳节点 xn 上旳近似值 yn，用以近似y(xn)
x0
y0
x1 f y(x), x dx
x0
x2 f y(x), x dx
x1
y(x1) f y(x1), x1 h
第17页
同样，在[x0,xn+1] ，积分采用矩形近似，得：
y(xn1) y0
f xn1
x0
y(x), x dx
y(xn ) f y(xn ), xn h
yn y(xn )
第5页
2、欧拉近似办法
2.1 简朴欧拉(L.Euler, 1707-1783)办法。
dy
dx
f
(y, x)
y(x0 ) y0
欧拉数值算法就是由初值通过递推求解，递推求解
就是从初值开始，后一种函数值由前一种函数值得到。核 心是构造递推公式。
y0 y1 y2 yn
第6页
i 1,2,...
第36页
没有一种算法可以有效地解决所有旳 ODE 问题，因此 MATLAB 提供了多种ODE函数。
函数 ODE类
特点
阐明
型
ode45
非刚性 单步法；4，5 阶 R-K 措施；合计 大部分场合旳首选措施
截断误差为 (△x)3
ode23
非刚性 单步法；2，3 阶 R-K 措施；合计 使用于精度较低旳情形
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用欧拉法求解常微分方程在数学和科学领域中，常微分方程是一种非常有用的工具，用于描述许多自然和物理现象。
MATLAB是一种强大的数学软件，可以用来解决许多数学问题。
本文将介绍如何使用欧拉法在MATLAB中求解常微分方程。
欧拉法是一种基本的数值方法，用于近似解决微积分方程问题。
该方法使用离散时间步长，将微积分方程转换成差分方程，并不断迭代求解。
欧拉法的实现非常简单，因此它很适合初学者。
下面是使用欧拉法在MATLAB中求解常微分方程的步骤：1. 定义常微分方程以 y' = -0.5y + 3sin(t) 为例，我们先定义常微分方程。
在MATLAB中，可以使用 anonymous functions 实现：dydt = @(t,y) -0.5*y + 3*sin(t);2. 定义时间范围和时间步长我们需要定义时间范围和时间步长，以便在一定时间范围内求解差分方程。
在这个例子中，我们定义时间范围为 0 到 10，并定义时间步长为 0.1：tspan = [0 10];h = 0.1;3. 定义初始条件我们需要定义初始条件，即 y(0) 的值。
在这个例子中，我们假设 y(0) = 1：y0 = 1;4. 求解差分方程现在我们可以使用欧拉法求解差分方程了。
在MATLAB中，可以使用 odeEuler 函数（需要自己编写）：[t,y] = odeEuler(dydt,tspan,y0,h);5. 可视化结果最后，我们可以将结果可视化，以便更好地理解求解过程。
在这个例子中，我们可以用 plot 函数将求解结果绘制出来：plot(t,y)xlabel('Time')ylabel('y(t)')title('Solution of y'' = -0.5y + 3sin(t) using Euler''s method')以上就是使用欧拉法在MATLAB中求解常微分方程的基本步骤。
matlab 数学应用 微分方程 常微分方程 非负ode解
matlab 数学应用微分方程常微分方程非负ode解
在MATLAB中解决微分方程，特别是常微分方程（ODEs），通常使用内置的ode45函数。
这个函数可以解决非负的常微分方程。
以下是一个简单的示例，说明如何使用ode45来解决一个简单的非负常微分方程：
假设我们要解决以下方程：
dy/dt = y - y^2
这个方程描述了一个在生物学或经济学中常见的模型，其中y表示某种数量，t表示时间。
以下是MATLAB代码：
y) y - y^2;
% 初始条件
y0 = 0.5; % 初始值
tspan = [0, 10]; % 时间范围
% 使用ode45求解
[t, y] = ode45(dydt, tspan, y0);
% 绘制结果
plot(t, y(:,1));
xlabel('Time');
ylabel('y(t)');
title('Solution of the ODE');
代码首先定义了微分方程（使用匿名函数@(t, y)）。
然后，它设置了初始条件和时间范围。
最后，它使用ode45来求解微分方程，并使用plot函数来绘制结果。
注意：ode45是默认使用四阶龙格-库塔法和五阶龙格-库塔法进行数值求解的，这两种方法都是相当稳定和可靠的。
但是，对于某些问题，可能需要尝试其他的数值方法或调整参数。
matlab用欧拉法求常微分方程初值
matlab用欧拉法求常微分方程初值用欧拉法求解常微分方程是一种常用的数值解法。
在数学和工程领域中，常微分方程是一类描述自然现象和物理过程的重要方程。
在实际问题中，我们往往难以得到准确的解析解，因此需要借助数值方法来近似求解。
欧拉法是其中一种简单而有效的数值解法。
让我们来了解一下常微分方程的基本概念。
常微分方程是指未知函数与其导数之间的关系式。
通常形式为dy/dx=f(x,y)，其中f(x,y)为已知的函数。
常微分方程的解就是满足该关系式的函数y(x)。
接下来，我们来看一下欧拉法的基本原理。
欧拉法的基本思想是将微分方程转化为差分方程，通过迭代计算来逼近解析解。
具体而言，我们将自变量x离散化为一系列的点，然后根据微分方程的导数定义，将微分项转化为差分项。
假设我们的求解区间为[x0,xn]，步长为h，那么我们可以得到近似解的递推公式为：y(i+1) = y(i) + h*f(x(i),y(i))其中，y(i)表示第i个点的函数值，x(i)表示第i个点的自变量值，f(x(i),y(i))表示在(x(i),y(i))处微分方程的导数值。
通过递推计算，我们可以得到离散点上的函数近似解。
当步长h足够小的时候，欧拉法可以得到较为精确的结果。
然而，需要注意的是，欧拉法的精度受到步长的限制，当步长过大时，误差会较大。
现在，我们来通过一个具体的例子来说明欧拉法的应用。
假设我们要求解如下的常微分方程：dy/dx = x^2其中，初始条件为y(0) = 1，求解区间为[0,1]。
我们可以将该微分方程转化为差分方程，并使用欧拉法进行求解。
我们将求解区间离散化，假设步长h=0.1，则我们可以得到离散点x0=0，x1=0.1，x2=0.2，...，x10=1。
然后，根据欧拉法的递推公式，我们可以得到近似解的计算过程如下：y(1) = y(0) + h*f(x(0),y(0))= 1 + 0.1*(0^2)= 1y(2) = y(1) + h*f(x(1),y(1))= 1 + 0.1*(0.1^2)= 1.001y(3) = y(2) + h*f(x(2),y(2))= 1.001 + 0.1*(0.2^2)= 1.004...y(10) = y(9) + h*f(x(9),y(9))= y(9) + 0.1*(0.9^2)通过逐步计算，我们可以得到离散点上的近似解。
Matlab ode函数 微分方程的数值解
odefun是函数句柄,可以是函数文件名,匿名函数句柄或内联函数名
tspan是区间[t0 tf]或者一系列散点[t0,t1,...,tf]
y0是初始值向量《Simulink与信号处理》
T返回列向量的时间点
Y返回对应T的求解列向量
[T,Y] = ode45(odefun,tspan,y0,options)
options是求解参数设置,可以用odeset在计算前设定误差,输出参数,事件等
[T,Y,TE,YE,IE] =ode45(odefun,tspan,y0,options)
在设置了事件参数后的对应输出
TE事件发生时间
YE事件解决时间
IE The index i of the event functionthat vanishes.
sol =ode45(odefun,[t0 tf],y0...)
sol结构体输出结果
编辑本段
示例
求解一阶常微分方程
odefun=@(t,y) (y+3*t)/t^2; %定义函数
tspan=[1 4]; %求解区间
y0=-2; %初值
[t,y]=ode45(odefun,tspan,y0);plot(t,y) %作图
title('t^2y''=y+3t,y(1)=-2,1<t<4') legend('t^2y''=y+3t')
xlabel('t')
ylabel('y')=y+3*t','y(1)=-2')
常微分方程初值问题的MATLAB解法
1用Matlab 求常微分方程(ODE)的初值问题(IVP)本节考虑一阶常微分方程000(,) ()u f t u t t Tu t u '=<≤⎧⎨=⎩ (1.1)的数值求解问题，包括算法公式及编程问题。
对一阶常微分方程组的初值问题010111120221220200120()(,,,,)(,,,,)() (,,,,)()m mm m m m m u t u u f t u u u u f t u u u u t u t t T u f t u u u u t u ⎧='=⎧⎪⎪'==⎪⎪<≤⎨⎨⎪⎪⎪⎪'==⎩⎩ (1.2)可以通过引入列向量0,,u u f化成类似(1.1)的形式000(,) ()u f t u t t Tu t u ⎧'=<≤⎪⎨=⎪⎩(1.3)其中1101122202120012()()(,,,,)()()(,,,,)(),,(,)()()(,,,,)m m m m mm u t u t f t u u u u t u t f t u u u u t u f t u u t u t f t u u u ⎛⎫⎛⎫⎛⎫⎪ ⎪ ⎪ ⎪ ⎪ ⎪=== ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎝⎭⎝⎭⎝⎭ (1.4)另外一个高阶常微分方程的初值问题()(1)0(1)(1)000000(,,,,,) (),(),,()m m m m u f t u u u u t t Tu t u u t u u t u ---'''⎫=<≤⎪⎬''===⎪⎭(1.5)也可以通过变换(1)123,,,,m m u u u u u u u u -'''==== 化成维微分方程组：1223112(,,,,)m mmm u u u u u uu f t u u u -'=⎧⎪'=⎪⎪⎨⎪'=⎪'=⎪⎩(1.6)我们在设计算法时通常先对一维一阶常微分方程(1.1)进行，然后再将这个算法写成适合求解(1.3)的向量形式，并以向量形式来进行编程。
matlab求解微分方程
matlab求解微分方程
Matlab是一种广受欢迎的工程数学软件，它可以解决复
杂的数学问题，其中包括求解常微分方程（ODEs）。
ODEs
在几何、物理和工程上都有很多应用，包括解决振动问题、热传递问题、控制系统设计问题等。
Matlab有一个内置的ODE解算器，可以根据给定的微分
方程求解。
Matlab中的ODE解算器可以分为两类：一类是ODE
45，它使用四阶Runge-Kutta法；另一类是ODE15s，它
使用单步积分器，如Adams-Bashforth和Adams-Moulton方法。
当使用Matlab求解ODE时，首先必须编写一个函数来定
义微分方程右端的函数，然后调用ODE解算器来求解，可以
指定边界条件、初始值、解算器类型等参数。
Matlab的可视化工具可以用来可视化求解的解，可以方
便地查看求解的结果，从而更好地理解问题。
总之，Matlab是一款功能强大的工程数学软件，它可以
很方便地求解ODE，可以用于各种科学和工程应用，是一款
非常强大的软件。
matlab常微分方程初值问题的数值解法
yn+1
h = y n + [ f ( x n, y n) + f ( x n+1 , yn+1 )] 2
y ( x 2 ) ≈ y ( x 0 ) + 2 h f ( x1 , y ( x1 ))
y′( x1 ) ≈ y ( x 2 ) − y ( x0 ) 2h
y n +1 = y n −1 + 2h f (xn , y n ) n = 1, 2, ⋯ 改进欧拉法 /* modified Euler’s method */
改 进 E u le r 公 式
数值分析
数值分析
以上两组公式都使用函数f ( x , y )在某些点上的 值的线性组合来计算y( xn +1 )的近似值yn +1。 的值，为一阶方法。 Euler公式： 公式 Euler公式：每步计算一次f ( x , y )的值，为一阶方法。 改进Euler公式：需计算两次f ( x , y )的值，二阶方法。 的值，二阶方法。 改进Euler公式 Euler公式：
例 ：用尤 拉公 式和 改进 的尤拉 公式 解初 值问 题 2x ' y = y− y y(0) = 1. (0 < x < 1);
数值分析
数值分析
解 ： 取 步 长 h = 0 .1， 尤 拉 公 式 为 ： yn+1 2 xn ). = yn + h( yn − yn
matlab的ode函数
matlab的ode函数在MATLAB中，ODE函数（ordinary differential equation）用于求解常微分方程（ordinary differential equations，ODEs）的数值解。
ODE函数在MATLAB的“ode”命令下调用，有多种不同类型的ODE求解器可供选择。
ODE函数的语法如下：[t, y] = ode45(odefun, tspan, y0)其中- “ode45”是一种常用的ODE求解器，可以用于求解非刚性的一阶或高阶常微分方程。
- “odefun”是一个函数句柄，表示待求解的ODE，其形式为dy/dt= f(t, y)。
该函数需要接受两个输入参数：自变量t和因变量y，并输出对应的导数值。
- “tspan”是一个包含两个元素的向量，表示自变量t的范围。
-“y0”是一个包含初始条件的列向量，表示因变量y在自变量t的初始值。
ODE函数的输出包括两个变量：- “t”是一个列向量，表示自变量t的离散值。
这些值等距地分布在tspan范围内。
-“y”是一个矩阵，每一列代表因变量y在相应t值处的数值解。
除了ode45之外，MATLAB还提供了其他常用的ODE求解器，包括ode23、ode113、ode15s、ode23s等。
这些求解器根据不同的算法和精度要求进行了优化，可根据具体问题的特点选择适当的求解器。
在使用ODE函数时，需要定义一个ODE函数句柄，作为输入参数传递给ODE求解器。
此函数句柄需要编写对应的ODE方程，并确保正确地输入输出导数值。
以下是一个示例ODE函数的编写过程：function dydt = myODE(t, y)dydt = 2 * y; % 示例ODE：dy/dt = 2*y然后，可以调用ODE求解器来求解该ODE：这将返回自变量t的离散值和对应的因变量y的数值解。
通过使用ODE函数，用户可以方便地在MATLAB环境中求解常微分方程，并获取其数值解。
常微分方程数值求解 MATLAB 求解
dsolve 举例
例：[x,y]=dsolve('Dx+5*x=0', 'Dy-3*y=0', ...
'x(0)=1', 'y(0)=1', 't') sol = dsolve('Dx+5*x=0', 'Dy-3*y=0', ...
'x(0)=1', 'y(0)=1', 't') 这里返回的 sol 是一个 结构类型 的数据 sol.x % 查看解函数 x(t) sol.y % 查看解函数 y(t)
常微分方程数值求解
—— MATLAB 求解
Matlab 解初值问题函数
用 Maltab自带函数 解初值问题
求解析解：dsolve
求数值解：
ode45、ode23、 ode113、ode23t、ode15s、 ode23s、ode23tb
2
符号求解
符号求解
dsolve
3
dsolve 求解析解
如果省略自变量，则默认自变量为 t
dsolve('Dy=2*x','x'); % dy/dx = 2x
dsolve('Dy=2*x');
% dy/dt = 2x
若找不到解析解，则提出警告，并返回空解。
5
dsolve 的使用
使用符号方程
导数：diff，如 diff(y)，diff(y,2) 等号：== 必须声明应变量与自变量!
求解析解：dsolve
y=dsolve('eq1','eq2', ... ,'cond1','cond2', ... ,'v')
第12章 MATLAB常微分方程(组)数值求解方程与方程组的数值解
ode15s求解起来效率不太高的刚性问题。
ode23t可以用来求解微分代数方程。
ode23tb 刚性
低
当方程是刚性的，并且求解要求精度不高
时可以使用。
2020/6/19
© 吴鹏, MATLAB从零到进阶.
常微分方程数值解
【例12.4-1】的结果图：
方 法 1计 算 结 果 图 1
方 法 2计 算 结 果 图 1
方 法 3计 算 结 果 图 1
0.8
0.8
0.8
0.6
0.6
0.6
0.4
0.4
0.4
0.2
0.2
0.2
0
0
0
-0.2
-0.4 0
2020/6/19
y1(t)
-0.2
y2(t)
y3(t)
-0.4
在某段时间区间内的变化来看。非刚性问题变化相对缓 慢，而刚性问题在某段时间内会发生剧烈变化，即很短 的时间内，解的变化巨大。对于刚性问题不适合用 ode45来求解，如果硬要用ode45来求解的话，达到指定 精度所耗费的时间往往会非常长 。
2020/6/19
© 吴鹏, MATLAB从零到进阶.
二、 非刚性问题举例
这些函数可以求解非刚性问题，刚性问题，隐式
微分方程，微分代数方程等初值问题，也可以求解 延迟微分方程，以及边值问题等。
2020/6/19
© 吴鹏, MATLAB从零到进阶.
二、初值问题求解函数
常微分方程数值解
1. 提供的函数
ode23, ode45, ode113, ode15s, ode23s, ode23t, ode23tb，这些函数统一 的调用格式如下：
常微分方程的数值解的matlab命令实现方法
常微分方程的数值解的matlab命令实现方法常微分方程的数值解在 MATLAB 中可以通过 ode 函数或 dsolve 函数进行求解。
其中，ode 函数可以求解一阶常微分方程，而 dsolve 函数可以求解二阶及以上的常微分方程。
下面是具体的实现方法:1. 一阶常微分方程的求解对于一阶常微分方程，可以使用 ode 函数求解。
假设我们要求解的常微分方程为:dx/dt = f(x, t)可以使用以下命令进行求解:y0 = [a, 0]; % 初值条件tspan = [0, 20]; % 时间区间[t, y] = ode45(@(t, y) odefun(t, y, a), tspan, y0); % 求解其中，odefun 函数用于定义常微分方程的解，它是一个自定义函数，其形式可以为:dy/dt = f(t, y)其中，dy 是 y 的求导，f(t, y) 是常微分方程的系数矩阵。
在 MATLAB 中，可以使用 dy[] 函数来计算 y 的求导，例如:dy = dy[](t, y);最后，使用 ode45 函数求解常微分方程的解，其中 tspan 是时间区间，y0 是初值条件。
2. 二阶常微分方程的求解对于二阶常微分方程，可以使用 dsolve 函数求解。
假设我们要求解的二阶常微分方程为:d2y/dt2 + p(t)dyy/dt + q(t)dy/dt + r(t)y = 0可以使用以下命令进行求解:syms t pqr;y0 = [a1, a2, a3]; % 初值条件[t, y] = dsolve(@(t, y) dy0(t, y), t, y0); % 求解其中，dy0 函数用于定义二阶常微分方程的解，其形式可以为:d2y/dt2 + p(t)dyy/dt + q(t)dy/dt + r(t)y = 0其中，d2y/dt2 是 y 的二阶求导，其它项是 y 的求导。
在 MATLAB 中，可以使用 dy0[] 函数来计算 y 的二阶求导。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
4
2.微分方程转换
好， 上面我们把 Matlab 中的常微分方程(ODE)的解算器讲解的差不多了， 下面我们就具体开始介绍如何使用上面 的知识吧！ 现实总是残酷的，要得到就必须先付出，不可能所有的 ODE 一拿来就可以直接使用，因此，在使用 ODE 解算 器之前，我们需要做的第一步，也是最重要的一步就是将微分方程组化成 Matlab 可接受的标准形式！ 如果 ODEs 由一个或多个高阶微分方程给出，则我们应先将其变换成一阶显式常微分方程组！ 下面我们以两个高阶微分方程构成的 ODEs 为例介绍如何将之变换成一个一阶显式常微分方程组。
【解】真是万幸，该 ODEs 已经帮我们完成了 step1，我们只需要完成 step2 和 step3 了
(1)选择一组状态变量
x1 = x, x2 = x ', x3 = y, x4 = y '
(2)写出所有状态变量的一阶微分表达式
6
x1 ' = x2 x2 ' = 2 x4 + x1 − µ* ( x1 + µ ) / r13 − µ ( x1 − µ * ) / r23 x3 ' = x4 x4 ' = −2 x2 + x3 − µ * x3 / r13 − µ x3 / r23
在求解 ODE 时，我们还会用到 deval()函数，deval 的作用就是通过结构体 solution 计算 t 对应 x 值，和 polyval 之类的很相似 输入参数： sol：就是上次调用 ode**函数得道的结构体解 xint：需要计算的点，可以是标量或者向量，但是必须在 tspan 范围内 该函数的好处就是如果我想知道 t=t0 时的 y 值，不需要重新使用 ode 计算，而直接使用上次计算的得到 solution 就可以
注意到，最高阶次的微分式的表达式直接就是 step1 中的微分方程 好，到此为止，一阶显式常微分方程组，变化顺利结束，接下来就是 Matlab 编程了。请记住上面的变化很重要， Matla 中所有微分方程的求解都需要上面的变换
下面通过一个实例演示 ODEs 的转换和求解
已知Apoll卫星的运动轨迹( x, y )满足下面的微分方程组 µ *( x + µ ) µ (x − µ* ) − x '' = 2 y '+ x − r13 r23 * y '' = −2 x '+ y − µ y − µ y r13 r23 其中： r1 = ( x + µ )2 + y 2 , r2 = ( x − µ * )2 + y 2 , µ * = 1 − µ , µ = 1/ 82.45 x(0) = 1.2, x '(0) = 0, y (0) = 0, y '(0) = −1.04935751
非刚性
低到高
刚性(stiff) 刚性 中等刚性 刚性
低到中 低 低 低
相关优化选项 参数 relto abstol normcontrol outputfcn outputsel refine stats nonnegative events maxstep inialstep jacobian jpattern vectorized mass mstatedependence ode45 √ ode23 √ ode113 √ ode15s √ ode23s √ ode23t √ ode23tb √
Matlab 中常微分方程数值解法讲解©
1.ODE 解算器简介 ..............................................................................................................................................................3 2.微分方程转换 ...................................................................................................................................................................5 3.刚性/非刚性问题 ..............................................................................................................................................................8 4.隐式微分方程(IDE) ........................................................................................................................................................10 5.微分代数方程(DAE).......................................................................................................................................................15 6.延迟微分方程(DDE).......................................................................................................................................................18 7.边值问题(BVP) ...............................................................................................................................................................20
然而现实总是残酷的，有时方程偏偏是隐式的，没法写成上面的样子，不用担心 Matlab 早就为我们想到了，这 个留在后面的隐式微分方程数值求解中再详细讲解！
step2.为每一阶微分式选择状态变量，最高阶除外
x1 = x, x2 = x ', x3 = x '', L, xm = x ( m −1) xm+1 = y, xm+1 = y ', xm + 2 = y '',L , xm + n = y ( n−1)
step1.将微分方程的最高阶变量移到等式的左边，其他移到右边，并按阶次从低到高排列，假如说两个高阶微
分方程最后能够显式的表达成如下所示：
x ( m ) = f (t , x, x ', x '', L , x ( m −1) , y , y ', L , y ( n−1) ) (n) ( m −1) , y , y ', L , y ( n−1) ) y = g (t , x, x ', x '', L , x
√ √ √ √
× × × ×
√ √ √ ×
√ × × ×
注：√*表示只能应用于没有质量矩阵的问题 输入参数： odefun：微分方程的句柄，必须写成 Matlab 规范格式，这个具体在后面讲解 tspab=[t0 tf]或者[t0,t1,…tf]：微分变量的范围，两者都是根据 t0 和 tf 的值自动选择步长，后只是前者返回所有计 算点的微分值，而后者只返回指定的点的微分值，一定要注意对于后者 tspan 必须严格单调，还有就是两者数据 存储时使用的内存不同(明显前者多)，其它没有任何本质的区别 y0=[y(0),y'(0),y"(0)…]：微分方程初值，依次输入所有状态变量的初值 options：微分优化参数，是一个结构体，使用 odeset 可以设置具体参数，详细内容查看帮助 输出参数： T：时间列向量 Y：二维数组，第 i 列表示第 i 个状态变量的值， 行数与 T 一致
从上面的变换，我们注意到，ODEs 中所有因变量的最高阶次之和就是需要的状态变量的个数，最高阶的微分式 (比如上面的 x-(m)和 y(n))不需要给它一个状态变量
step3.根据上面选用的状态变量，写出所有状态变量的一阶微分的表达式
5
x '1 = x2 x '2 = x3 x '3 = x4 M x 'm = f (t , x1 , x2 , x3 ,L , xm + n ) x 'm+1 = xm + 2 M x 'm+ n = g (t , x1 , x2 , x3 ,L , xm+ n )