ODE常微分方程
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
end
设t0 = 0(年)， x0 = 1(万台)， x = 100 (万台) ， = 0.01(年-1 万台-1)， 可用下 列命令作出8年内电饭锅销量预测图形：
»fplot('100/(1+99*exp(-0.01*100*x))',[0,8]) 可见短期预报二个模型相近，但作为长 期预报，后者较前者合理。当然后者也 有不尽合理之处，比如x难以确定，未 考虑产品更新换代等。
end
例2 解方程y’ = y-2t/y, y(0)=1, 0<t<1 解 先写M函数eg5_2fun.m » [t,y]=euler('eg5_2fun',[0,1],1,0.1)
四、使用MATLAB命令 1、数值解 [tout,yout] = ode45('yprime', tspan, y0) 用法与euler相同。 若无输出参数，则作出图形。 ode23 与ode45类似只是精度低一些。
end
2、初等积分法 3、常系数线性微分方程 线性常微分方程的解为一个特解和相应 的齐次微分方程通解的叠加。 齐次微分方程的解可用特征根法求得 例1 求x’’+ 0.2 x’+3.92x = 0的通解 解 特征方程为2 + 0.2 +3.92=0 »roots([1 0.2 3.92 0] 求得共轭复根 +i=-0.11.9774i, 通解为 x(t) = Aet cos(t) +Bet sin(t)
Euler格式
y (t k 1 ) y (t k ) y '(t k ) h k=0,1,2…
end
yk 1 yk hf (t k , yk )
M函数euler.m给出Euler法计算程序 使用格式为 [tout,yout] = euler('ypfun', tspan, y0,h) ypfun: 表示f(t, y)的M文件名 tspan=[t0, tf]: 表示自变量初值t0和终值tf y0: 表示初值向量y0，h是步长。 输出列向量tout: 表示节点 (t0 , t1 , … , tn)' 输出矩阵yout: 表示数值解，每一列对应 y的一个分量
一、引例 导弹系统的改进 海军方面要求改进现有的舰对舰导弹系 统。目前的电子系统能迅速测出敌舰的种 类、位置以及敌舰行驶速度和方向，且导 弹自动制导系统能保证在发射后任一时刻 都能对准目标。根据情报，这种敌舰能在 我军舰发射导弹后T分钟作出反应并摧毁 导弹。现在要求改进电子导弹系统使能自 动计算出敌舰是否在有效打击范围之内。
end
在线算法：对于测定的d 和，可用
dx dt
dy dt
b dy 2 1 ( ) dx
b dx 2 1 ( ) dy
b at sin y (t ) 2 1 ( ) d at cos x(t )
b d at cos x(t ) 2 1 ( ) at sin y (t )
内，可以发射。
end
二、数学理论复习: 常微分方程 1、微分方程的概念
常微分方程: f (t,y,y’,y’’,…,y(n))=0 微分方程组: 联系一些未知函数的一组微 分方程 线性常微分方程: y(n) + a1 (t) y(n-1) + … + an-1 (t) y’ + an (t) y = b(t) 若ai (t) (i =1, …,n) 与t无关, 称为常系数的 若b(t)=0,称为齐次的
end
三、微分方程数值解法：Euler法 t0 t t f y ' (t ) f (t , y (t )) y (t 0 ) y0 数值解法:寻求解y(t)在一系列离散节点 t0 < t1 < …< tn <tf 上的近似值yk (k=0,1,…n)。 hk = tk+1 －tk 为步长，通常取为常量h 。 Euler法：在节点处用差商近似代替导数
end
五、实验例题 例5（引例）在导弹系统中设 a=90km/h, b=450km/h, T=0.1h. 求d, 的有效范围? 解 有两个极端情形容易算出。 若 =0, 即敌舰正好背向行驶, 即x轴正向。 那么导弹直线飞行, 击中时间 t=d/(b-a)<T 得d=T(b-a)=36km。 若 =, 即迎面驶来, 类似有d=T(a+b)=54km 一般地, 有36<d<54。
b d at cos x(t ) 2 1 ( ) at sin y (t )
end
初始条件为x(0)=0, y(0)=0, 对于给定的
a,b,d, 进行计算。当x(t)满足
x(t) + > d + a t cos
则认为已击中目标。
这里代表允许的误差，因为敌舰是有一
定大小的。如果t < T,则敌舰在打击范围
end
设x为全部需要量，那么销售速度与当时 的潜在需要量 (x - x) 成正比，则有方程：
dx x ( x x ) dt
其中为比例常数。可用dsolve » dsolve(' Dx=a*x*(x1-x)','x(t0)=x0') 解得 x
x (t )
x 1 ( 1) exp( x (t t0 )) x0
计算出t。如d=50, =/2，写出M函数 eg5_5fun.m 用euler即得 x=44.2893
end
由于在T小时内，横坐标没有突破x=50, 所以敌舰不在有效打击范围，应等近一些 再发射。 离线算法：首先对于所有可能的d和，计 算击中所需时间，从而对不同 ，得d的临 界值。具体应用时直接查表判断。 x(t) + > d + a t cos 取 =0.1，编写m脚本文件eg5_5.m 运行得临界曲线。使用时查询即可。
end
例6 经调查发现，电饭锅销售速度与当时 的销量成正比。现在我们来建立一个数学 模型以预测销量。设x(t)表示t时刻的销量， x0为初始时刻t0的销量，那么有方程
dx kx dt 其中k为常数。解得 x(t) = x0 e k(t-t0)。 当k > 0, t时，x(t)， 这对于销售初 期可认为是合适的，长期显然不合适。
end
d
end
为了保持对准目标，导弹轨迹切线方向 应为 dy at sin y (t ) dx d at cos x(t ) 由上面两个方程得下列微分方程
dx dt
dy dt
b Βιβλιοθήκη Baiduy 2 1 ( ) dx
b dx 2 1 ( ) dy
b at sin y (t ) 2 1 ( ) d at cos x(t )
end
设我舰发射导弹时位置在坐标原点，敌 舰在x轴正向d（km）处, 其行驶速度为 a (km/h), 方向与x轴夹角为 , 导弹飞行线速 度b(km/h) 。设t 时刻时导弹位置为 (x(t)，y(t)) , 那么
dx 2 dy 2 ( ) ( ) b dt dt
易知 t 时刻敌舰位置 为 (d+atcos,atsin)。
3、刚性方程组解法 刚性方程组解法ode15s使用格式同ode45
y1 ' 0.01y1 99.99 y2 例4 100 y2 y2 ' y1 (0) 2 y2 (0) 1
解 先将方程写为M函数eg5_4fun.m »[t,y]=ode15s('eg5_4fun',[0,400],[2,1]’); »plot(t,y);
end
2、符号微分方程解析解 s=dsolve(‘方程1’,‘方程2’,,‘初始条 件 1’,‘初始条件2’,,‘自变量’) 均用字符串方式表示,自变量缺省值为t, 导数用D表示,2阶导数用D2表示,以此类推。 s返回解析解, 方程组情形, s为符号结构。
例3 (1)求y'=ay+b的通解；(2)求解例2 (3)高阶方程 y''=cos(2x)-y, y(0)=1, y'(0)=0 (4)方程组 f '=f+g, g'=-f+g, f(0)=1, g(0)=2end
设t0 = 0(年)， x0 = 1(万台)， x = 100 (万台) ， = 0.01(年-1 万台-1)， 可用下 列命令作出8年内电饭锅销量预测图形：
»fplot('100/(1+99*exp(-0.01*100*x))',[0,8]) 可见短期预报二个模型相近，但作为长 期预报，后者较前者合理。当然后者也 有不尽合理之处，比如x难以确定，未 考虑产品更新换代等。
end
例2 解方程y’ = y-2t/y, y(0)=1, 0<t<1 解 先写M函数eg5_2fun.m » [t,y]=euler('eg5_2fun',[0,1],1,0.1)
四、使用MATLAB命令 1、数值解 [tout,yout] = ode45('yprime', tspan, y0) 用法与euler相同。 若无输出参数，则作出图形。 ode23 与ode45类似只是精度低一些。
end
2、初等积分法 3、常系数线性微分方程 线性常微分方程的解为一个特解和相应 的齐次微分方程通解的叠加。 齐次微分方程的解可用特征根法求得 例1 求x’’+ 0.2 x’+3.92x = 0的通解 解 特征方程为2 + 0.2 +3.92=0 »roots([1 0.2 3.92 0] 求得共轭复根 +i=-0.11.9774i, 通解为 x(t) = Aet cos(t) +Bet sin(t)
Euler格式
y (t k 1 ) y (t k ) y '(t k ) h k=0,1,2…
end
yk 1 yk hf (t k , yk )
M函数euler.m给出Euler法计算程序 使用格式为 [tout,yout] = euler('ypfun', tspan, y0,h) ypfun: 表示f(t, y)的M文件名 tspan=[t0, tf]: 表示自变量初值t0和终值tf y0: 表示初值向量y0，h是步长。 输出列向量tout: 表示节点 (t0 , t1 , … , tn)' 输出矩阵yout: 表示数值解，每一列对应 y的一个分量
一、引例 导弹系统的改进 海军方面要求改进现有的舰对舰导弹系 统。目前的电子系统能迅速测出敌舰的种 类、位置以及敌舰行驶速度和方向，且导 弹自动制导系统能保证在发射后任一时刻 都能对准目标。根据情报，这种敌舰能在 我军舰发射导弹后T分钟作出反应并摧毁 导弹。现在要求改进电子导弹系统使能自 动计算出敌舰是否在有效打击范围之内。
end
在线算法：对于测定的d 和，可用
dx dt
dy dt
b dy 2 1 ( ) dx
b dx 2 1 ( ) dy
b at sin y (t ) 2 1 ( ) d at cos x(t )
b d at cos x(t ) 2 1 ( ) at sin y (t )
内，可以发射。
end
二、数学理论复习: 常微分方程 1、微分方程的概念
常微分方程: f (t,y,y’,y’’,…,y(n))=0 微分方程组: 联系一些未知函数的一组微 分方程 线性常微分方程: y(n) + a1 (t) y(n-1) + … + an-1 (t) y’ + an (t) y = b(t) 若ai (t) (i =1, …,n) 与t无关, 称为常系数的 若b(t)=0,称为齐次的
end
三、微分方程数值解法：Euler法 t0 t t f y ' (t ) f (t , y (t )) y (t 0 ) y0 数值解法:寻求解y(t)在一系列离散节点 t0 < t1 < …< tn <tf 上的近似值yk (k=0,1,…n)。 hk = tk+1 －tk 为步长，通常取为常量h 。 Euler法：在节点处用差商近似代替导数
end
五、实验例题 例5（引例）在导弹系统中设 a=90km/h, b=450km/h, T=0.1h. 求d, 的有效范围? 解 有两个极端情形容易算出。 若 =0, 即敌舰正好背向行驶, 即x轴正向。 那么导弹直线飞行, 击中时间 t=d/(b-a)<T 得d=T(b-a)=36km。 若 =, 即迎面驶来, 类似有d=T(a+b)=54km 一般地, 有36<d<54。
b d at cos x(t ) 2 1 ( ) at sin y (t )
end
初始条件为x(0)=0, y(0)=0, 对于给定的
a,b,d, 进行计算。当x(t)满足
x(t) + > d + a t cos
则认为已击中目标。
这里代表允许的误差，因为敌舰是有一
定大小的。如果t < T,则敌舰在打击范围
end
设x为全部需要量，那么销售速度与当时 的潜在需要量 (x - x) 成正比，则有方程：
dx x ( x x ) dt
其中为比例常数。可用dsolve » dsolve(' Dx=a*x*(x1-x)','x(t0)=x0') 解得 x
x (t )
x 1 ( 1) exp( x (t t0 )) x0
计算出t。如d=50, =/2，写出M函数 eg5_5fun.m 用euler即得 x=44.2893
end
由于在T小时内，横坐标没有突破x=50, 所以敌舰不在有效打击范围，应等近一些 再发射。 离线算法：首先对于所有可能的d和，计 算击中所需时间，从而对不同 ，得d的临 界值。具体应用时直接查表判断。 x(t) + > d + a t cos 取 =0.1，编写m脚本文件eg5_5.m 运行得临界曲线。使用时查询即可。
end
例6 经调查发现，电饭锅销售速度与当时 的销量成正比。现在我们来建立一个数学 模型以预测销量。设x(t)表示t时刻的销量， x0为初始时刻t0的销量，那么有方程
dx kx dt 其中k为常数。解得 x(t) = x0 e k(t-t0)。 当k > 0, t时，x(t)， 这对于销售初 期可认为是合适的，长期显然不合适。
end
d
end
为了保持对准目标，导弹轨迹切线方向 应为 dy at sin y (t ) dx d at cos x(t ) 由上面两个方程得下列微分方程
dx dt
dy dt
b Βιβλιοθήκη Baiduy 2 1 ( ) dx
b dx 2 1 ( ) dy
b at sin y (t ) 2 1 ( ) d at cos x(t )
end
设我舰发射导弹时位置在坐标原点，敌 舰在x轴正向d（km）处, 其行驶速度为 a (km/h), 方向与x轴夹角为 , 导弹飞行线速 度b(km/h) 。设t 时刻时导弹位置为 (x(t)，y(t)) , 那么
dx 2 dy 2 ( ) ( ) b dt dt
易知 t 时刻敌舰位置 为 (d+atcos,atsin)。
3、刚性方程组解法 刚性方程组解法ode15s使用格式同ode45
y1 ' 0.01y1 99.99 y2 例4 100 y2 y2 ' y1 (0) 2 y2 (0) 1
解 先将方程写为M函数eg5_4fun.m »[t,y]=ode15s('eg5_4fun',[0,400],[2,1]’); »plot(t,y);
end
2、符号微分方程解析解 s=dsolve(‘方程1’,‘方程2’,,‘初始条 件 1’,‘初始条件2’,,‘自变量’) 均用字符串方式表示,自变量缺省值为t, 导数用D表示,2阶导数用D2表示,以此类推。 s返回解析解, 方程组情形, s为符号结构。
例3 (1)求y'=ay+b的通解；(2)求解例2 (3)高阶方程 y''=cos(2x)-y, y(0)=1, y'(0)=0 (4)方程组 f '=f+g, g'=-f+g, f(0)=1, g(0)=2end