道客巴巴微分方程

数 学 建 模 与 数 学 实 验后勤工程学院数学教研室 微 分 方 程 实 验 目 的实 验 内 容 2、 学 会 用 Matlab求 微 分 方 程 的 数 值 解 .1、 学 会 用 Matlab求 简 单 微 分 方 程 的 解 析 解 .1、 求 简 单 微 分 方 程 的 解 析 解 .4、 实 验 作 业 .2、 求 微 分 方 程 的 数 值 解 .3、 数 学 建 模 实 例 求 微 分 方 程 的 数 值 解（ 一 ） 常 微 分 方 程 数 值 解 的 定 义（ 二 ） 建 立 数 值 解 法 的 一 些 途 径（ 三 ） 用 Matlab软 件 求 常 微 分 方 程 的 数 值 解返 回 1、 目 标 跟 踪 问 题 一 ： 导 弹 追 踪 问 题 2、 目 标 跟 踪 问 题 二 ： 慢 跑 者 与 狗3、 地 中 海 鲨 鱼 问 题 返 回数 学 建 模 实 例 微 分 方 程 的 解 析 解 求 微 分 方 程 （ 组 ） 的 解 析 解 命 令 :dsolve(方 程 1, 方 程 2,方 程 n, 初 始 条 件 , 自 变 量 ) 记 号 : 在 表 达 微 分 方 程 时 ， 用 字 母 D 表 示 求 微 分 ， D2、 D3 等 表 示 求 高 阶 微 分 .任 何 D 后 所 跟 的 字 母 为 因 变 量 ， 自 变 量 可 以 指 定 或 由 系 统 规 则 选 定 为 确 省 . 例 如 ， 微 分 方 程 02 2 dx yd 应 表 达 为 ： D2y=0. 例 1 求 21 udt du 的 通 解 . 解 输 入 命 令 ： dsolve(Du=1+u2,t) To Matlab（ ff1） 结 果 ： u = tg(t-c) 例 2 求 微 分 方 程 的 特 解 . 15)0(,0)0( 029422 yy ydxdydxyd 解 输 入 命 令 : y=dsolve(D2y+4*Dy+29*y=0,y(0)=0,Dy(0)=15,x)结 果 为 : y =3e-2xsin（ 5x） To Matlab（ ff2） 例 3 求 微 分 方 程 组 的 通 解 . zyxdtdz zyxdtdy zyxdtdx 244 354 332解 输 入 命 令 ： x,y,z=dsolve(Dx=2*x-3*y+3*z,Dy=4*x-5*y+3*z,Dz=4*x-4*y+2*z, t)； x=simple(x) % 将 x化 简 y=simple(y) z=simple(z)结 果 为 ： x = (c 1-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 To Matlab（ ff3） 返 回 微 分 方 程 的 数 值 解（ 一 ） 常 微 分 方 程 数 值 解 的 定 义 在 生 产 和 科 研 中 所 处 理 的 微 分 方 程 往 往 很 复 杂 且 大 多得 不 出 一 般 解 。 而 在 实 际 上 对 初 值 问 题 ， 一 般 是 要 求 得到 解 在 若 干 个 点 上 满 足 规 定 精 确 度 的 近 似 值 ， 或 者 得 到一 个 满 足 精 确 度 要 求 的 便 于 计 算 的 表 达 式 。因 此 ， 研 究 常 微 分 方 程 的 数 值 解 法 是 十 分 必 要 的 。 。的 相 应 近 似 值 求 出 准 确 值，值 处 ， 即 对的 若 干 离 散 的 开 始其 数 值 解 是 指 由 初 始 点，：对 常 微 分 方 程 nn nyyxyxy xxxxx y ,y )(,),( ),y(x x )y(x y)f(x,y 212 1210 000 返 回 （ 二 ） 建 立 数 值 解 法 的 一 些 途 径 001i )y(x y)f(x,y ,1,2,1,0 , x ynihxi 解 微 分 方 程 ：可 用 以 下 离 散 化 方 法 求设 1、 用 差 商 代 替 导 数 若 步 长 h较 小 ， 则 有 h xyhxyxy )()()( 故 有 公 式 ： 1-n,0,1,2,i )( ),(00 1 xyy yxhfyy iiii此 即 欧 拉 法 。 2、 使 用 数 值 积 分对 方 程 y=f(x,y), 两 边 由 xi到 xi+1积 分 ， 并 利 用 梯 形 公 式 ， 有 ：)(,()(,(2)(,()()( 1111 1 iiiiiixxii xyxfxyxfxxdttytfxyxy ii实 际 应 用 时 ， 与 欧 拉 公 式 结 合 使 用 ： ,2,1,0 ),(),(2 ),( )( 11)1( 1)0( 1 kyxfyxfhyy yxhfyy kiiiiiki iiii 的 计 算 。然 后 继 续 下 一 步 ，取时 ，当 满 足，对 于 已 给 的 精 确 度 ）（ y y 2i 111i)( 1)1( 1 kikiki yyy 此 即 改 进 的 欧 拉 法 。故 有 公 式 ： )( ),(),(200 111 xyy yxfyxfhyy iiiiii 3、 使 用 泰 勒 公 式 以 此 方 法 为 基 础 ， 有 龙 格 -库 塔 法 、 线 性 多 步 法 等 方法 。4、 数 值 公 式 的 精 度 当 一 个 数 值 公 式 的 截 断 误 差 可 表 示 为 O（ hk+1） 时（ k为 正 整 数 ， h为 步 长 ） ， 称 它 是 一 个 k阶 公 式 。k越 大 ， 则 数 值 公 式 的 精 度 越 高 。欧 拉 法 是 一 阶 公 式 ， 改 进 的 欧 拉 法 是 二 阶 公 式 。龙 格 -库 塔 法 有 二 阶 公 式 和 四 阶 公 式 。线 性 多 步 法 有 四 阶 阿 达 姆 斯 外 插 公 式 和 内 插 公 式 。返 回 （ 三 ） 用 Matlab软 件 求 常 微 分 方 程 的 数 值 解t， x=solver（ f,ts,x0,options）ode45 ode23 ode113ode15sode23s 由 待 解方 程 写成 的 m-文 件 名 ts=t0， tf，t0、 tf为 自变 量 的 初值 和 终 值 函 数 的初 值ode23： 组 合 的 2/3阶 龙 格 -库 塔 -芬 尔 格 算 法ode45： 运 用 组 合 的 4/5阶 龙 格 -库 塔 -芬 尔 格 算 法自 变量 值 函 数值 用 于 设 定 误 差 限 (缺 省 时 设 定 相 对 误 差 10 -3, 绝 对 误 差 10-6),命 令 为 ： options=odeset（ reltol,rt,abstol,at） , rt， at： 分 别 为 设 定 的 相 对 误 差 和 绝 对 误 差 . 1、 在 解 n个 未 知 函 数 的 方 程 组 时 ， x0和 x均 为 n维 向 量 ，m-文 件 中 的 待 解 方 程 组 应 以 x的 分 量 形 式 写 成 . 2、 使 用Matlab软 件 求 数 值 解 时 ， 高 阶 微 分 方 程 必 须等 价 地 变 换 成 一 阶 微 分 方 程 组.注 意 : 例 4 0)0(;2)0( 0)1(1000 222 xx xdtdxxdtxd解: 令 y1=x，y2=y1 则 微 分 方 程 变 为 一 阶 微 分 方 程 组 ： 0)0(,2)0( )1(1000 21 12212 21 yy yyyy yy1、 建 立m-文 件vdp1000.m如 下 ： function dy=vdp1000(t,y) dy=zeros(2,1); dy(1)=y(2); dy(2)=1000*(1-y(1)2)*y(2)-y(1); 2、 取 t 0=0， tf=3000， 输 入 命 令 ： T,Y=ode15s(vdp1000,0 3000,2 0); plot(T,Y(:,1),-)3、 结 果 如 图 0 500 1000 1500 2000 2500 3000-2.5-2-1.5-1-0.500.511.52 To Matlab（ ff4） 例 5 解 微 分 方 程 组 . 1)0(,1)0(,0)0( 51.0 321 213 312 321 yyy yyy yyy yyy解 1、 建 立m-文 件rigid.m如 下 ： function dy=rigid(t,y) dy=zeros(3,1); dy(1)=y(2)*y(3); dy(2)=-y(1)*y(3); dy(3)=-0.51*y(1)*y(2);2、 取 t 0=0， tf=12， 输 入 命 令 ： T,Y=ode45(rigid,0 12,0 1 1); plot(T,Y(:,1),-,T,Y(:,2),*,T,Y(:,3),+)3、 结 果 如 图 To Matlab（ ff5）0 2 4 6 8 10 12-1-0.8-0.6-0.4-0.2 0 0.2 0.4 0.6 0.8 1 图 中 ， y1的 图 形 为 实 线 ， y2的 图 形 为 “ *”线 ， y3的 图 形 为 “ +”线 .返 回 导 弹 追 踪 问 题 设 位 于 坐 标 原 点 的 甲 舰 向 位 于 x轴 上 点 A(1, 0)处 的 乙 舰发 射 导 弹 ， 导 弹 头 始 终 对 准 乙 舰 .如 果 乙 舰 以 最 大 的 速 度v0(是 常 数 )沿 平 行 于 y轴 的 直 线 行 驶 ， 导 弹 的 速 度 是 5v0， 求导 弹 运 行 的 曲 线 方 程 .又 乙 舰 行 驶 多 远 时 ， 导 弹 将 它 击 中 ？解 法 一 （ 解 析 法 ） 假 设 导 弹 在 t 时 刻 的 位 置 为 P(x(t), y(t)， 乙 舰 位 于 ),1( 0tvQ . 由 于 导 弹 头 始 终 对 准 乙 舰 ， 故 此 时 直 线 PQ 就 是 导 弹 的 轨 迹 曲 线 弧 OP 在 点 P 处 的 切 线 ， 即 有 xytvy 1 0 即 yyxtv )1(0 (1) 又 根 据 题 意 ， 弧 OP 的 长 度 为 AQ 的 5 倍 ， 即 tvdxyx 00 2 51 (2) 由 (1),(2)消 去 t整 理 得 模 型 : (3) 151)1( 2yyx 初 值 条 件 为 : 0)0( y 0)0( y 解 即 为 导 弹 的 运 行 轨 迹 : 245)1(125)1(85 5 654 xxy 当 1x 时 245y ， 即 当 乙 舰 航 行 到 点 )245 ,1( 处 时 被 导 弹 击 中 . 被 击 中 时 间 为 : 00 245vvyt . 若 v0=1, 则 在 t=0.21 处 被 击 中 .To Matlab(chase1)轨 迹 图 见 程 序 chase1 解 法 二 (数 值 解 )1.建 立m-文 件eq1.m function dy=eq1(x,y) dy=zeros(2,1); dy(1)=y(2); dy(2)=1/5*sqrt(1+y(1)2)/(1-x); 2. 取x 0=0，xf=0.9999， 建 立 主 程 序ff6.m如 下 ： x0=0，xf=0.9999 x,y=ode15s(eq1,x0 xf,0 0); plot(x,y(:,1),b.) hold on y=0:0.01:2; plot(1,y,b*) 结 论 : 导 弹 大 致 在 （1，0.2） 处 击 中 乙 舰 To Matlab(ff6)2151)1( yyx )1/(151 212 21 xyy yy 令 y1=y,y2=y1， 将 方 程 （ 3） 化 为 一 阶 微 分 方 程 组 。 解 法 三 (建 立 参 数 方 程 求 数 值 解 ) 设 时 刻 t乙 舰 的 坐 标 为 (X(t)， Y(t)， 导 弹 的 坐 标 为 (x(t)， y(t). 1 设 导 弹 速 度 恒 为 w， 则 222 )()( wdtdydtdx （ 1） 2. 由 于 弹 头 始 终 对 准 乙 舰 ， 故 导 弹 的 速 度 平 行 于 乙 舰 与 导 弹 头 位 置 的 差 向 量 ， 即 ： yY xXdtdydtdx ， 0 （ 2） 消 去 得 ： )()()( )()()( 22 22 yYyYxX wdtdy xXyYxX wdtdx （ 3）3 因 乙 舰 以 速 度 v0沿 直 线 x=1运 动 ， 设 v0=1， 则 w=5， X=1， Y=t 因 此 导 弹 运 动 轨 迹 的 参 数 方 程 为 ： 0)0(,0)0( )()()1( 5 )1()()1( 5 22 22yx ytytxdtdy xytxdtdx4. 解 导 弹 运 动 轨 迹 的 参 数 方 程建 立 m-文 件 eq2.m如 下 ： function dy=eq2(t,y) dy=zeros(2,1); dy(1)=5*(1-y(1)/sqrt(1-y(1)2+(t-y(2)2); dy(2)=5*(t-y(2)/sqrt(1-y(1)2+(t-y(2)2); 取t 0=0，tf=2， 建 立 主 程 序chase2.m如 下 ： t,y=ode45(eq2,0 2,0 0); Y=0:0.01:2; plot(1,Y,-)， hold on plot(y(:,1),y(:,2),*)To Matlab(chase2) 5. 结 果 见 图 1导 弹 大 致 在 （1，0.2） 处 击 中 乙 舰 ， 与 前 面 的 结 论 一 致.图 1 图 2 返 回 在chase2.m中 ， 按 二 分 法 逐 步 修 改tf， 即 分 别 取tf=1，0.5,0.25,直 到tf=0.21时 ， 得 图2.结 论 ： 时 刻 t=0.21时 ， 导 弹 在 （ 1， 0.21） 处 击 中 乙 舰 。To Matlab(chase2) 慢 跑 者 与 狗 一 个 慢 跑 者 在 平 面 上 沿 椭 圆 以 恒 定 的 速 率 v=1跑 步 ,设 椭圆 方 程 为 : x=10+20cost, y=20+5sint. 突 然 有 一 只 狗 攻 击 他 . 这 只 狗从 原 点 出 发 ,以 恒 定 速 率 w跑 向 慢 跑 者 ,狗 的 运 动 方 向 始 终 指 向 慢 跑 者 .分 别 求 出 w=20,w=5时 狗 的 运 动 轨 迹 .1. 模 型 建 立设 时 刻 t慢 跑 者 的 坐 标 为 (X(t)， Y(t)， 狗 的 坐 标 为 (x(t)， y(t). 则 X=10+20cost, Y=20+15sint, 狗 从 (0,0)出 发 ,与 导 弹 追 踪 问 题 类 似 ， 建 立 狗的 运 动 轨 迹 的 参 数 方 程 : 0)0( ,0)0( )sin1520()sin1520()cos2010( )cos2010()sin1520()cos2010( 22 22yx ytytxt wdtdy xtytxt wdtdx 2. 模 型 求 解(1) w=20时 ,建 立 m-文 件 eq3.m如 下 : function dy=eq3(t,y) dy=zeros(2,1); dy(1)=20*(10+20*cos(t)-y(1)/sqrt (10+20*cos(t)-y(1)2+(20+15*sin(t)-y(2)2); dy(2)=20*(20+15*sin(t)-y(2)/sqrt (10+20*cos(t)-y(1)2+(20+15*sin(t)-y(2)2);取t0=0，tf=10， 建 立 主 程 序chase3.m如 下 ： 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),*) 在chase3.m， 不 断 修 改 tf的 值 ,分 别 取 tf=5, 2.5, 3.5,至 3.15时 ,狗 刚 好 追 上 慢 跑 者 . To Matlab(chase3) 建 立 m-文 件 eq4.m如 下 : function dy=eq4(t,y) dy=zeros(2,1); dy(1)=5*(10+20*cos(t)-y(1)/sqrt (10+20*cos(t)-y(1)2+(20+15*sin(t)-y(2)2); dy(2)=5*(20+15*sin(t)-y(2)/sqrt (10+20*cos(t)- y(1)2+(20+15*sin(t)-y(2)2);取t0=0，tf=10， 建 立 主 程 序chase4.m如 下 ： t0=0;tf=10; t,y=ode45(eq4,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),*) 在chase3.m， 不 断 修 改 tf的 值 ,分 别 取 tf=20, 40, 80,可 以 看 出 ,狗 永 远 追 不 上 慢 跑 者 . To Matlab(chase4) (2) w=5时 返 回 地 中 海 鲨 鱼 问 题 意 大 利 生 物 学 家 Ancona曾 致 力 于 鱼 类 种 群 相 互 制 约 关系 的 研 究 ， 他 从 第 一 次 世 界 大 战 期 间 ,地 中 海 各 港 口 捕 获 的几 种 鱼 类 捕 获 量 百 分 比 的 资 料 中 ， 发 现 鲨 鱼 等 的 比 例 有 明显 增 加 （ 见 下 表 ） ， 而 供 其 捕 食 的 食 用 鱼 的 百 分 比 却 明 显下 降 .显 然 战 争 使 捕 鱼 量 下 降 ， 食 用 鱼 增 加 ， 鲨 鱼 等 也 随 之增 加 ， 但 为 何 鲨 鱼 的 比 例 大 幅 增 加 呢 ？ 他 无 法 解 释 这 个 现 象 ， 于 是 求 助 于 著 名 的 意 大 利 数 学家 V.Volterra， 希 望 建 立 一 个 食 饵 捕 食 系 统 的 数 学 模 型 ，定 量 地 回 答 这 个 问 题 . 年 代 1914 1915 1916 1917 1918 百 分 比 11.9 21.4 22.1 21.2 36.4 年 代 1919 1920 1921 1922 1923 百 分 比 27.3 16.0 15.9 14.8 19.7 1 符 号 说 明 ：)(1 tx 食 饵 在 t 时 刻 的 数 量 ； )(2 tx 捕 食 者 在 t 时 刻 的 数 量 ；1r 食 饵 独 立 生 存 时 的 增 长 率 ； 2r 捕 食 者 独 自 存 在 时 的 死 亡 率 ；1 捕 食 者 掠 取 食 饵 的 能 力 ； 2 食 饵 对 捕 食 者 的 供 养 能 力 . e捕 获 能 力 系 数 2 基 本 假 设 ： （ 1） 食 饵 由 于 捕 食 者 的 存 在 使 增 长 率 降 低 ， 假 设 降 低 的 程 度 与 捕 食 者 数 量 成 正 比 ； （ 2） 捕 食 者 由 于 食 饵 为 它 提 供 食 物 的 作 用 使 其 死 亡 率 降 低 或 使 之 增 长 ， 假 定 增 长 的 程 度 与 食 饵 数 量 成 正 比 。 3 模 型 建 立 与 求 解 模 型 （ 一 ） 不 考 虑 人 工 捕 获)( 21111 xrxdtdx )( 12222 xrxdtdx 该 模 型 反 映 了 在 没 有 人 工 捕 获 的 自 然 环 境 中 食 饵 与 捕 食 者 之 间 的 制 约关 系 ， 没 有 考 虑 食 饵 和 捕 食 者 自 身 的 阻 滞 作 用 ， 是 Volterra提 出 的 最 简 单 的模 型 . 针 对 一 组 具 体 的 数 据 用 Matlab 软 件 进 行 计 算 . 设 食 饵 和 捕 食 者 的 初 始 数 量 分 别 为 101 )0( xx ， 202 )0( xx 对 于 数 据 2,25,02.0,5.0,1.0,1 20102211 xxrr ，t 的 终 值 经 试 验 后 确 定 为 15， 即 模 型 为 ： 2)0(,25)0( )02.05.0( )1.01( 21 122 211 xx xxx xxx首 先 ， 建 立 m-文 件shier.m如 下 ： function dx=shier(t,x) dx=zeros(2,1); dx(1)=x(1)*(1-0.1*x(2); dx(2)=x(2)*(-0.5+0.02*x(1);其 次 ， 建 立 主 程 序 shark.m如 下 ： t,x=ode45(shier,0 15,25 2); plot(t,x(:,1),-,t,x(:,2),*) plot(x(:,1),x(:,2) To Matlab(shark) 相 图 ),( 21 xx 为 ： 0 5 10 150 10 20 30 40 50 60 70 80 90 100 0 20 40 60 80 1000 5 10 15 20 25 30 数 值 解 如 下 图 ： )(1 tx 为 实 线 ， )(2 tx 为 “ *” 线 . 求 解 结 果 ： 左 图 反 映 了 x 1（ t） 与 x2（ t） 的 关 系 。 可 以 猜 测 ： x1（ t） 与 x2（ t） 都 是 周 期 函 数 。 模 型 （ 二 ） 考 虑 人 工 捕 获 设 表 示 捕 获 能 力 的 系 数 为 e， 相 当 于 食 饵 的 自 然 增 长 率由 r1 降 为 r1-e， 捕 食 者 的 死 亡 率 由 r2 增 为 r2+e )( )( 12222 21111 xerxdtdx xerxdtdx 20,250,02.0,5.0,1.0,1 212211 ）（）（仍 取 xxrr 设 战 前 捕 获 能 力 系 数 e=0.3, 战 争 中 降 为 e=0.1, 则 战 前 与 战 争 中 的 模 型 分 别 为 : 2)0(,25)0( )02.08.0( )1.07.0( 21 122 211 xx xxdtdx xxdtdx 2)0(,25)0( )02.06.0( )1.09.0( 21 122 211 xx xxdtdx xxdtdx 模 型 求 解 :1、 分 别 用 m-文 件 shier1.m和 shier2.m定 义 上 述 两 个 方 程2、 建 立 主 程 序 shark1.m, 求 解 两 个 方 程 ， 并 画 出 两 种 情 况 下鲨 鱼 数 在 鱼 类 总 数 中 所 占 比 例 x2(t)/x1(t)+x2(t)To Matlab(shark1) 0 5 10 150 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 实 线 为 战 前 的 鲨鱼 比 例 ， “ *”线 为战 争 中 的 鲨 鱼 比 例结 论 ： 战 争 中 鲨 鱼 的 比 例 比 战 前 高 ！ 返 回 实 验 作 业 1. 一 个 小 孩 借 助 长 度 为 a的 硬 棒 ， 拉 或 推 某 玩 具 .此 小 孩 沿 某 曲线 行 走 ， 计 算 并 画 出 玩 具 的 轨 迹 . 2. 讨 论 资 金 积 累 、 国 民 收 入 、 与 人 口 增 长 的 关 系 .（ 1） 若 国 民 平 均 收 入 x与 按 人 口 平 均 资 金 积 累 y成 正 比 ， 说 明仅 当 总 资 金 积 累 的 相 对 增 长 率 k大 于 人 口 的 相 对 增 长 率 r时 ， 国民 平 均 收 入 才 是 增 长 的 .（ 2） 作 出 k(x)和 r(x)的 示 意 图 ， 分 析 人 口 激 增 会 引 起 什 么 后 果 . 返 回