matlab常微分方程课件

1.483240
0.8 1.649783
1.612452
1.0 1.784770
1.732051
误差
0.008602 0.016572 0.025726 0.037331 0.052719
与准确解 y 12x相比，可看出欧拉公式的计算结
果精度很差.
上页 下页
欧拉公式具有明显的几何意义, 就是用折线近似代
或表示为下列平均化形式
yp yc
yn yn
hf hf
( x n , y n ), ( x n1 , y p ),
f(x ,y 1 ) f(x ,y 2 ) L y 1 y 2 .
它表明f满足利普希茨条件，且L的大小反映了右端函 数f关于y变化的快慢，刻画了初值问(1.1)式和(1.2)式 是否为好条件. 这在数值求解中也是很重要的.
上页 下页
虽然求解常微分方程有各种各样的解析方法，但 解析方法只能用来求解一些特殊类型的方程，实际问 题中归结出来的微分方程主要靠数值解法.
利普希茨常数(简称Lips.常数).
上页 下页
定理1 设f在区域D={(x, y)|axb,yR}上连续, 关 于y满足利普希茨条件，则对任意x0[a, b], y0R，常 微分方程初值问题(1.1)式和(1.2)式当x[a, b]时存在 唯一的连续可微解y(x) .
解的存在唯一性定理是常微分方程理论的基本内 容，也是数值方法的出发点，此外还要考虑方程的解 对扰动的敏感性，它有以下结论.
可得欧拉法(2.1)的公式误差为
y (x n 1 ) y n 1 h 2 2y (n ) h 2 2y (x n ), (2 .3 )
称为此方法的局部截断误差.
上页 下页

取t0=0，tf=10，建立程序如下： t0=0;tf=10; [t,y]=ode45('eq3',[t0 tf],[0 0]); T=0:0.1:2*pi; X=10+20*cos(T); Y=20+15*sin(T); plot(X,Y,'-') hold on plot(y(:,1),y(:,2),'*')
注意: 1、在解n个未知函数的方程组时，x0和x均为n维向量，
m-文件中的待解方程组应以x的分量形式写成.
2、使用Matlab软件求数值解时，高阶微分方程必须 变换成等价的一阶微分方程组.
例4
d2x dt 2
1000(1
x2
)
dx dt
x
0
x(0) 2; x'(0) 0
解: 令 y1=x，y2=y1’
由待解 方程写 成的m-
ts=[t0，tf]， t0、tf为自
函数的 初值
变量的初
文件名 值和终值
ode23：组合的2/3阶龙格-库塔-芬尔格算法 ode45：运用组合的4/5阶龙格-库塔-芬尔格算法
用于设定误差限(缺省时设定相对误差10-3, 绝对误差10-6), 命令为：options=odeset（’reltol’,rt,’abstol’,at）, rt，at：分别为设定的相对误差和绝对误差.
z = (-C1+exp(4*t)*C1-2*exp(4*t)+C2+exp(4*t)*C3)*exp(-2*t)
用Matlab软件求常微分方程的数值解
[t，x]=solver（’f’,ts,x0,options）
自变 函数 量值 值
ode45 ode23 ode113 ode15s ode23s

科学计算与MATLAB
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

odeplot 随计算过程，显示解向
量
属性名 OutputSel
Refine
2020/4/10
参数设置
取值
含义
若不使用缺省设置，则
有效值： 正整数向 量 缺省值：[]
OutputFcn所表现的是那 些正整数指定的解向量 中的分量的曲线或数据 。若为缺省值时，则缺 省地按上面情形进行操
作
有效值： 正整数k>1 缺省值：k
y1 (0 ) y 0
y0
y
2
(
0
)
y1
yn
(
0
)
yn
• （3）根据（1）与（2）的结果，编写能计 算
• 导数的M-函数文件odefile。
• （4）将文件odefile与初始条件传递给求解
• 器Solver中的一个，运行后就可得到ODE
• 的、在指定时间区间上的解列向量y（其中包 2含020/4y/10及不同阶的导数）。
参数设置
属性名 NormControl Events
取值
含义
有效值： on、off 缺省值：
off
为‘on’时，控制解向量 范数的相对误差，使每 步计算中，满足： norm(e)<=max(RelTol*n orm(y),AbsTol)
有效值： 为‘on’时，返回相应的 on、off 事件记录
2020/4/10
例3
x' x 2
x(0)
1
创建函数function2, 保存在function2.m中
function f=function2(t,x)
f=-x.^2;
在命令窗口中执行 >> [t,x]=ode45('function2',[0,1],1); >> plot(t,x,'-',t,x,'o'); >> xlabel('time t0=0,tt=1'); >> ylabel('x values x(0)=1');

0
, y(t),r 均可测出。
可算出白铅中铅的衰变率 y 0 ，再于当时的矿物 比较，以鉴别真伪。 矿石中铀的最大含量可能 2~3%，若白铅中铅210 每分钟衰变超过3 万个原子，则矿石中含铀量超 过 4%。
测定结果与分析
画名 Emmaus的信徒们 洗足 钋210衰变原子数 镭226衰变原子数
更少量的镭(Ra226)。白铅是由铅金属产生的，而
铅金属是经过熔炼从铅矿中提取来出的。当白铅
从处于放射性平衡状态的矿中提取出来时， Pb210 的绝大多数来源被切断，因而要迅速蜕变，直到 Pb210与少量的镭再度处于放射平衡，这时Pb210 的蜕变正好等于镭蜕变所补足的为止。
铀238
T 45 亿年
4
6
Байду номын сангаас
5 5 x 1 y ( 1 , ) 当 时 ， 即 当 乙 舰 航 行 到 点 处 时 被 导 弹 击 中 . 24 24 y5 t . 被 击 中 时 间 为 : 若 v = 1 , 则 在 t = 0 . 2 1 处 被 击 中 . 0 v v 0 24 0
二 范. 梅格伦（Van Meegren） 伪造名画案
( t t ) 0
N t t0 ln 0 N 1
半衰期 T
1

ln 2
, N(t) 能测出或算出，只要知道 N 0 就可算出
断代。 这正是问题的难处，下面是间接确定N 0 的方法。
油画中的放射性物质 白铅（铅的氧化物）是油画中的颜料之一，应
用已有2000余年，白铅中含有少量的铅(Pb210)和
镭226
T 1600 年
（无放射性）
铅206
钋210

作线性组合得到平均斜率，由此得到更高阶的精
度，这就是龙格-库塔方法的基本思路。
在MATLAB软件中含有数值求解的系统函数，其
实现原理就是龙格－库塔方法。
其中参数k >0，m=18。
2020/11/12
8
2.2.2 利用平衡与增长式
许多研究对象在数量上常常表现出某种不变的特性， 如封闭区域内的能量、货币量等。利用变量间的平衡 与增长特性，可分析和建立有关变量间的相互关系。
此类建模方法的关键是分析并正确描述基本模型的右 端，使平衡式成立。
例2.2 战斗模型：两方军队交战，希望为这场战斗建 立一个数学模型，应用这个模型达到如下目的：
“速率”、“增长”、“衰变”、“边际的”等常涉及到
导数。
我们熟悉的速度公式：dy v 就是一个简单的一阶微分方
程。
dt
微分方程是指含有导数或微分的等式。 一般形式： F(x,y,y,,y(n))0
或y： (n) f(x,y,y,,y(n1)).
常用的建立微分方程的方法有：运用已知物理定律；利用 平衡与增长式；运用微元法；应用分析法。
y2(k+1)=(y2(k)+h*x1(k+1)+h)/(1+h);
y3(k+1)=(y3(k)+(h/2)*(-y3(k)+x1(k)+x1(k+1)+2))/(1+h/2);
end
x=0:0.1:1;
y=x+exp(-x);
x1=x1(1:11),y=y(1:11),y1=y1(1:11),y2=y2(1:11),y3=y3( 1:11),
202092040捕食者的存在使食饵的增长率降低假设x1降低的程度与捕食者数量x2成正比即食饵对捕食者的数量x2起到增长的作用其程度与食饵数量x1成正比即dtdxdtdxdtdxdtdx生产理论把企业仅仅抽象为一个生产函数一种投入产出关系一个追求利润最大化的黑匣子它没有讨论企业内部是如何配置资源的企业是如何组织生产的企业和市场的关系如何各自的边界在哪里

π
y′′ − µ 1 − y 2 y′ + y = 0 y (0) = 1, y′(0 ) = 0, µ = 2
函数ode23和ode45是对一阶常微分方程组设计的，因此，对高阶常微 和 是对一阶常微分方程组设计的， 函数 是对一阶常微分方程组设计的 因此， 分方程，需先将它转化为一阶常微分方程组，即状态方程。 分方程，需先将它转化为一阶常微分方程组，即状态方程。 则可写出Van der Pol 方程的状态方程形式： 方程的状态方程形式： 令 x1 = y , x2 = y ′, 则可写出
求解器 Ode23 Ode45 Ode113 Ode23t Ode15s Ode23s Ode23tb ode15i 采用方法 2-3阶Runge-Kutta算法，低精度 算法， 阶 算法 4-5阶Runge-Kutta算法，中精度 阶 算法， 算法 Adams算法，精度可到 10 −3 ~ 10 −6 算法， 算法 梯形算法 Gear’s方向数值微分算法，中精度 方向数值微分算法， 方向数值微分算法 2阶Rosebrock算法，低精度 阶 算法， 算法 梯形算法，低精度 梯形算法， 可变秩方法 适用场合 非刚性 非刚性 非刚性， 非刚性， 适度刚性 刚性 刚性 刚性 完全隐式微分方程
∫
1
0
ln xdx .
（1）建立被积函数文件 ）建立被积函数文件feln.m. function y=feln(x) y=exp(x).*log(x); （2）调用数值积分函数 求定积分。 ）调用数值积分函数quadgk求定积分。 求定积分 format long; I=quadgk(@ feln,0,1) 3.梯形积分法 梯形积分法 在MATLAB中，对由表格形式定义的函数关系的求定积分问题用梯 中 形积分函数trapz.该函数调用格式如下。 该函数调用格式如下。 形积分函数 该函数调用格式如下 ●T=trapz(Y).例如：trapz([1:5;2:6]’) 例如： 例如 ●T=trapz(X,Y).