matlab常微分方程 ppt课件
合集下载
MATLAB第十三讲常微分方程初值问题数值解法 ppt课件
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 )
称为此方法的局部截断误差.
上页 下页
matlab常微分方程数值解法精品PPT课件
科学计算与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
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常微分方程 ppt课件
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');
第8章MATALAB微分方程-PPT精选文档
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教学第二章 微分方程PPT课件
作线性组合得到平均斜率，由此得到更高阶的精
度，这就是龙格-库塔方法的基本思路。
在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生产理论把企业仅仅抽象为一个生产函数一种投入产出关系一个追求利润最大化的黑匣子它没有讨论企业内部是如何配置资源的企业是如何组织生产的企业和市场的关系如何各自的边界在哪里
MATLAB课件 第6讲
π
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).
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
2020/4/10
求解具体ODE的基本过程：
• （2）运用数学中的变量替换：
• yn=y(n-1),yn-1=y(n-2),…,y2=y‘,y1=y，
• 把高阶（大于2阶）的方程（组）写成一阶
•
微分y 方程组yyy：n12
f1 ( t , y )
f2 (t,
y
)
fn
(t,
y
)
2020/4/10
2020/4/10
二、常微分方程（组）数值解
• Matlab专门用于求解常微分方程的函 • 数，主要采用Runge-Kutta方法： • ode23, ode45, ode113, ode15s, • ode23s, ode23t, ode23tb
2020/4/10
二、常微分方程（组）数值解
• [T,Y] = solver(odefun,tspan,y0) • [T,Y] = solver(odefun,tspan,y0,options) • [T,Y] = • solver(odefun,tspan,y0,options,p1,p2…)
对于标量t与列向量y，函数f=odefun(t,y) 必须返回一f(t,y)的列向量f。
解矩阵Y中的每一行对应于返回的时间列向 量T中的一个时间点。
要获得问题在其他指定时间点t0,t1,t2,… 上的解，则令tspan=[t0,t1,t2,…,tf]（要求 是单调的）。
2020/4/10
[T,Y] = solver(odefun,tspan,y0,options)
2020/4/10
参数说明：
• （4）y0 包含初始条件的向量。 • （5）options 用命令odeset设置的可选 • 积分参数. • （6）p1,p2,… 传递给函数odefun的可选 • 参数。
2020/4/10
[T,Y] = solver(odefun,tspan,y0)
在区间tspan=[t0,tf]上，从t0到tf，用 初始条件y0求解显式微分方程y’=f(t,y)。
y1 (0 ) y 0
y0
y
2
(
0
)
y1
yn
(
0
)
yn
• （3）根据（1）与（2）的结果，编写能计 算
• 导数的M-函数文件odefile。
• （4）将文件odefile与初始条件传递给求解
将参数p1,p2,p3,..等传递给函数 odefun，再进行计算。若没有参数设置，则 令options=[]。
2020/4/10
求解具体ODE的基本过程：
• （1）根据问题所属学科中的规律、定律、 • 公式，用微分方程与初始条件进行描述。 • F(y,y’,y’’,…,y(n),t) = 0 • y(0)=y0,y’(0)=y1,…,y(n-1)(0)=yn-1 • 而y=[y;y(1);y(2);…,y(m-1)]， • n与m可以不等
表示 y(x) xa b
y(x) xc d
y(x) xe f
2020/4/10
• (4)若边界条件少于方程（组）的阶数 ，则返回的结果r中会出现任意常数C1
• ，(5)Cd2s，olv…e命；令最多可以接受12个输入参量
（包括方程组与定解条件个数，当然我们 可以做到输入的方程个数多于12个，只要 将多个方程置于一字符串内即可）。
• ans=
• cos(a*x)
• >> dsolve(‘D2y = -a^2*y’, ‘y(0) = 1,Dy(pi/a) = 0’)
2020/4/10
u v例2
v
u
• >> [u,v] = dsolve('Du=v,Dv=u')
• u= • C1*exp(-t)+C2*exp(t) • V= • -C1*exp(-t)+C2*exp(t)
用参数options（用命令o来自eset生成）设置 属性（代替了缺省的积分参数），再进行操作 。常用的属性包括相对误差值RelTol（缺省值 为1e-3）与绝对误差向量AbsTol（缺省值为每 一元素为1e-6）。
2020/4/10
[T,Y] = solver(odefun,tspan,y0,options,p1,p2…)
• 科学技术和工程中许多问题是用微分方程的形 • 式建立数学模型，因此微分方程的求解有很实 • 际的意义。 • 一、常微分方程（组）的符号解 • 二、常微分方程（组）数值解
2020/4/10
一、常微分方程（组）的符号解
• 函数 dsolve 格式： • r = dsolve('eq1,eq2,…’,'cond1,cond2,…','v')
2020/4/10
参数说明：
• （1）solver为命令 • Ode45,ode23,ode113,ode15s, • ode23s,ode23t,ode23tb之一。 • （2）odefun 为常微分方程y’=f(x,y)， • 或为包含一混合矩阵的方程(x,y)*y’=f(x,y). • （3）tspan 积分区间（即求解区间）的向 • 量tspan=[t0,tf]。要获得问题在其他指定 • 时间点t0,t1,t2,…上的解，则令 • tspan=[t0,t1,t2,…,tf] • （要求是单调的）。
2020/4/10
• (6)若没有给定输出参量，则在命令窗口显 示解列表。若该命令找不到解析解，则返 回一警告信息，同时返回一空的sym对象。 这时，用户可以用命令ode23或ode45求解 方程组的数值解。
2020/4/10
例1
y a 2 y
y(0) 1
y (
/
a)
0
• >> dsolve('D2y = -a^2*y,y(0) = 1,Dy(pi/a) = 0','x')
说明 (1)对给定的常微分方程（组）eq1,eq2,… 中指定的符号自变量v，与给定的边界条件 和初始条件cond1,cond2,….求符号解（即解 析解）r； (2)若没有指定变量v，则缺省变量为t；
2020/4/10
(3)在微分方程（组）的表达式eq中，大写 字母D表示对自变量(设为x)的微分算子： D=d/dx，D2=d2/dx2，…。微分算子D后面的 字母则表示为因变量，即待求解的未知函 数。初始和边界条件由字符串表示： y(a)=b，Dy(c)=d，D2y(e)=f，等等，分别
求解具体ODE的基本过程：
• （2）运用数学中的变量替换：
• yn=y(n-1),yn-1=y(n-2),…,y2=y‘,y1=y，
• 把高阶（大于2阶）的方程（组）写成一阶
•
微分y 方程组yyy：n12
f1 ( t , y )
f2 (t,
y
)
fn
(t,
y
)
2020/4/10
2020/4/10
二、常微分方程（组）数值解
• Matlab专门用于求解常微分方程的函 • 数，主要采用Runge-Kutta方法： • ode23, ode45, ode113, ode15s, • ode23s, ode23t, ode23tb
2020/4/10
二、常微分方程（组）数值解
• [T,Y] = solver(odefun,tspan,y0) • [T,Y] = solver(odefun,tspan,y0,options) • [T,Y] = • solver(odefun,tspan,y0,options,p1,p2…)
对于标量t与列向量y，函数f=odefun(t,y) 必须返回一f(t,y)的列向量f。
解矩阵Y中的每一行对应于返回的时间列向 量T中的一个时间点。
要获得问题在其他指定时间点t0,t1,t2,… 上的解，则令tspan=[t0,t1,t2,…,tf]（要求 是单调的）。
2020/4/10
[T,Y] = solver(odefun,tspan,y0,options)
2020/4/10
参数说明：
• （4）y0 包含初始条件的向量。 • （5）options 用命令odeset设置的可选 • 积分参数. • （6）p1,p2,… 传递给函数odefun的可选 • 参数。
2020/4/10
[T,Y] = solver(odefun,tspan,y0)
在区间tspan=[t0,tf]上，从t0到tf，用 初始条件y0求解显式微分方程y’=f(t,y)。
y1 (0 ) y 0
y0
y
2
(
0
)
y1
yn
(
0
)
yn
• （3）根据（1）与（2）的结果，编写能计 算
• 导数的M-函数文件odefile。
• （4）将文件odefile与初始条件传递给求解
将参数p1,p2,p3,..等传递给函数 odefun，再进行计算。若没有参数设置，则 令options=[]。
2020/4/10
求解具体ODE的基本过程：
• （1）根据问题所属学科中的规律、定律、 • 公式，用微分方程与初始条件进行描述。 • F(y,y’,y’’,…,y(n),t) = 0 • y(0)=y0,y’(0)=y1,…,y(n-1)(0)=yn-1 • 而y=[y;y(1);y(2);…,y(m-1)]， • n与m可以不等
表示 y(x) xa b
y(x) xc d
y(x) xe f
2020/4/10
• (4)若边界条件少于方程（组）的阶数 ，则返回的结果r中会出现任意常数C1
• ，(5)Cd2s，olv…e命；令最多可以接受12个输入参量
（包括方程组与定解条件个数，当然我们 可以做到输入的方程个数多于12个，只要 将多个方程置于一字符串内即可）。
• ans=
• cos(a*x)
• >> dsolve(‘D2y = -a^2*y’, ‘y(0) = 1,Dy(pi/a) = 0’)
2020/4/10
u v例2
v
u
• >> [u,v] = dsolve('Du=v,Dv=u')
• u= • C1*exp(-t)+C2*exp(t) • V= • -C1*exp(-t)+C2*exp(t)
用参数options（用命令o来自eset生成）设置 属性（代替了缺省的积分参数），再进行操作 。常用的属性包括相对误差值RelTol（缺省值 为1e-3）与绝对误差向量AbsTol（缺省值为每 一元素为1e-6）。
2020/4/10
[T,Y] = solver(odefun,tspan,y0,options,p1,p2…)
• 科学技术和工程中许多问题是用微分方程的形 • 式建立数学模型，因此微分方程的求解有很实 • 际的意义。 • 一、常微分方程（组）的符号解 • 二、常微分方程（组）数值解
2020/4/10
一、常微分方程（组）的符号解
• 函数 dsolve 格式： • r = dsolve('eq1,eq2,…’,'cond1,cond2,…','v')
2020/4/10
参数说明：
• （1）solver为命令 • Ode45,ode23,ode113,ode15s, • ode23s,ode23t,ode23tb之一。 • （2）odefun 为常微分方程y’=f(x,y)， • 或为包含一混合矩阵的方程(x,y)*y’=f(x,y). • （3）tspan 积分区间（即求解区间）的向 • 量tspan=[t0,tf]。要获得问题在其他指定 • 时间点t0,t1,t2,…上的解，则令 • tspan=[t0,t1,t2,…,tf] • （要求是单调的）。
2020/4/10
• (6)若没有给定输出参量，则在命令窗口显 示解列表。若该命令找不到解析解，则返 回一警告信息，同时返回一空的sym对象。 这时，用户可以用命令ode23或ode45求解 方程组的数值解。
2020/4/10
例1
y a 2 y
y(0) 1
y (
/
a)
0
• >> dsolve('D2y = -a^2*y,y(0) = 1,Dy(pi/a) = 0','x')
说明 (1)对给定的常微分方程（组）eq1,eq2,… 中指定的符号自变量v，与给定的边界条件 和初始条件cond1,cond2,….求符号解（即解 析解）r； (2)若没有指定变量v，则缺省变量为t；
2020/4/10
(3)在微分方程（组）的表达式eq中，大写 字母D表示对自变量(设为x)的微分算子： D=d/dx，D2=d2/dx2，…。微分算子D后面的 字母则表示为因变量，即待求解的未知函 数。初始和边界条件由字符串表示： y(a)=b，Dy(c)=d，D2y(e)=f，等等，分别