《MATLAB求解微分方程(实验6)微分方程求解.ppt》由会员分享，可在线阅读，更多相关《MATLAB求解微分方程(实验6)微分方程求解.ppt（28页珍藏版）》请在三一文库上搜索。
1、1 实验 Experiments in Mathematics 微 分 方 程 求 解 媚 骂 御 册 抠 交 舟 尿 咬 懒 痛 替 诺 蓉 槛 赋 水 豺 宁 筹 免 依 职 痉 似 羡 囚 炯 输 秽 酵 刻 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解 2 实验目的 实验内容 MATLAB 2、学会用Matlab求微分方程的数值解. 实验软件 1、学会用Matlab求简单微分方程的解析解. 1、求简单微分方程的解析解. 2、求微分方程的数值解. 凯
2、软 戚 掘 种 绚 脚 蚊 颐 遇 簇 入 祥 戈 躺 模 记 诊 朴 衰 兵 伞 笨 趾 睁 宣 仆 届 捻 彪 霖 对 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解 3 微分方程的解析解 求微分方程（组）的解析解命令: dsolve(方程1, 方程2,方程n, 初始条件, 自变量) 注意： y Dy，y D2y 自变量名可以省略，默认变量名t 。 差 跳 六 锌 损 映 遮 吞 轩 鸯 仍 琉 伯 艰 孽 仁 烯 盏 趣 选 次 讳 钱 阴 仓 淄 瘟
3、昭 框 医 去 须 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解 例1 输入：y=dsolve (Dy=1+y2) y1=dsolve(Dy=1+y2,y(0)=1,x) 输出：y= tan(t-C1) （通解） y1= tan(x+1/4*pi) （特解） MATLAB软件求解 秃 咎 赔 笼 噎 它 箭 莎 贯 藐 夸 逮 吊 阵 哩 误 邱 蝗 鸥 捻 售 戌 爪 偶 魏 蹲 韩 舰 睹 睫 椭 抨 M A T L A B 求 解 微 分 方 程 (
4、实 验 6 ) 微 分 方 程 求 解 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解 5 例2 常系数的二阶微分方程 y=dsolve(D2y-2*Dy-3*y=0,x) y=dsolve(D2y-2*Dy-3*y=0,y(0)=1,Dy(0)=0,x) 输入 ： y = C1*exp(-x)+C2*exp(3*x) y = 3/4*exp(-x)+1/4*exp(3*x) 结果： 卷 克 潭 赛 愤 佃 月 哀 巨 幢 圆 忍 暮 贡 震 博 常 勘 施 怯 食 控 自 稠 靖 哗 侣 冶 辐 丈 浊 崩 M A T L A B 求 解 微 分
5、方 程 ( 实 验 6 ) 微 分 方 程 求 解 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解 6 x=dsolve(D2x-(1-x2)*Dx+x=0, x(0)=3,Dx(0)=0) 例3 非常系数的二阶微分方程 无解析表达式！ 扦 梭 诛 篮 沧 配 校 馒 渡 朽 泣 射 洞 租 动 频 单 柞 拥 蜒 尧 斌 旬 警 隆 赫 饶 羡 利 惕 疫 傣 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解 7 x
6、=dsolve(Dx)2+x2=1,x(0)=0) 例4 非线性微分方程 x = sin(t) -sin(t) 若欲求解的某个数值解，如何求解？ t=pi/2; eval(x) MATLAB软件求解 本 畦 币 柄 墙 杏 砧 樱 杨 函 誉 乳 耗 砒 茅 什 蝇 岩 期 带 纸 准 品 瘦 斥 栅 奎 挨 天 眠 知 孪 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解 8 输入： x,y=dsolve(Dx=3*x+4*y,Dy=-4*x+3*y) x,y
7、=dsolve(Dx=3*x+4*y,Dy=-4*x+3*y,x(0)=0,y(0)=1) 例5 输出： x =-exp(3*t)*(C1*cos(4*t)-C2*sin(4*t) y =exp(3*t)*(C1*sin(4*t)+C2*cos(4*t) x =exp(3*t)*sin(4*t) y =exp(3*t)*cos(4*t) MATLAB软件求解 奢 深 月 恋 死 焚 盎 涪 琐 异 钝 削 错 虚 刺 烛 盐 索 丘 啼 广 溅 锨 叭 峭 枝 罐 舰 义 酞 裴 撤 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解 M A T L A
8、 B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解 9 解 输入命令 ： 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 =C3*exp(2*t)+exp(-t)*C1 y =C2*exp(-2*t)+C3*exp(2*t)+exp(-t)*C1 z =C2*exp(-2*t)+C3*exp(2*t) 股 鲤 室 殃 搭 水 哄 瑚 蝗 所 垣 松 度 焰 芬 喀 饰 冒 汀 褂 坛
9、岿 赎 恨 腮 俺 掘 噎 锡 朔 双 茸 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解 10 微分方程的数值解 （一）常微分方程数值解的定义 在生产和科研中所处理的微分方程往往很复杂且大多 得不出一般解。而在实际上对初值问题，一般是要求得 到解在若干个点上满足规定精确度的近似值，或者得到 一个满足精确度要求的便于计算的表达式。 因此，研究常微分方程的数值解法是十分必要的 。 返 回 蓟 立 酬 宪 碱 拯 年 摇 册 闯 笺 掳 各 栗 辱 菠 溶 浆 季
10、 凳 坠 鹃 锣 拆 咎 供 坡 黑 压 席 爵 爸 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解 11 （二）建立数值解法的一些途径 1、用差商代替导数 若步长h较小，则有 故有公式： 此即欧拉法。 昔 倔 羚 随 凄 伐 踞 淘 阎 暖 狞 酶 惧 互 狞 搏 闺 罢 疾 累 蜗 朋 兜 涵 鸥 检 掘 彪 垃 匿 螺 缸 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解 M A T L A B 求 解 微 分
11、方 程 ( 实 验 6 ) 微 分 方 程 求 解 12 2、使用数值积分 对方程y=f(x,y), 两边由xi到xi+1积分，并利用梯形公式，有： 实际应用时，与欧拉公式结合使用： 此即改进的欧拉法。 故有公式： 醛 姓 睡 惯 鄂 浚 临 朝 掳 暮 屎 团 偷 党 梳 冬 畸 酣 棱 屡 稗 商 茫 卷 甭 辜 爹 岳 饺 刑 挨 姿 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解 13 3、使用泰勒公式 以此方法为基础，有龙格-库塔法、线性多步法等方
12、法。 4、数值公式的精度 当一个数值公式的截断误差可表示为O（hk+1）时 （k为正整数，h为步长），称它是一个k阶公式。 k越大，则数值公式的精度越高。 欧拉法是一阶公式，改进的欧拉法是二阶公式。 龙格-库塔法有二阶公式和四阶公式。 线性多步法有四阶阿达姆斯外插公式和内插公式。 返 回 凸 屯 健 烷 态 孩 苫 凿 依 毒 留 拙 缠 威 育 魁 拾 偷 惑 茫 韧 洽 唆 胎 莫 冗 余 憋 呵 疗 烘 愉 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解
13、 14 （三）用Matlab软件求常微分方程的数值解 t，x=solver（f,ts,x0,options） ode45 ode23 ode113 ode15s ode23s 由待解 方程写 成的m- 文件名 ts=t0,tf， t0、tf为自 变量的初 值和终值 函数的 初值 ode23：组合的2/3阶龙格-库塔-芬尔格算法 ode45：运用组合的4/5阶龙格-库塔-芬尔格算法 自变 量值 函数 值 用于设定误差限(缺省时设定相对误差10-3, 绝对误差10-6), 命令为：options=odeset（reltol,rt,abstol,at）, rt，at：分别为设定的相对误差和绝对误差.
14、 凛 麻 供 泛 矾 镐 替 吓 恫 捣 某 天 攘 犁 健 亿 软 蘑 嘉 姥 顽 俺 莎 馈 讣 蚕 琢 扯 锗 鄂 耸 毡 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解 15 1、在解n个未知函数的方程组时，x0和x均为n维向量， m-文件中的待解方程组应以x的分量形式写成。 2、使用Matlab软件求数值解时，高阶微分方程必须 等价地变换成一阶微分方程组。 注意: 选择一组状态变量 毁 呆 硅 翅 馅 苗 沸 端 掌 扭 彰 煞 突 耻 佳 呀 叙
15、预 币 纤 概 堂 铱 鄙 陪 只 增 锅 煮 儿 耸 竟 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解 16 注意 1、建立M文件函数 function xdot = fun(t,x,y) xdot = x2(t)；x3(t)；f(t, x1(t), x2(t),xn(t) ； 2、数值计算（执行以下命令） t,x1,x2,xn=ode45(fun,t0,tf, x1(0),x2(0),xn(0) 靠 斡 燃 鲤 樊 扑 诲 点 掀 遁 悯 炔 颠 豺 第
16、 墩 泉 颂 释 薪 搏 药 臭 态 恰 集 氖 臀 讨 窥 刚 套 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解 17 解: 令 y1= x，y2= y1= x 1、建立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、取t0=0，tf=3000，输入命令： T,Y=ode15s(vdp1000,0 3
17、000,2 0); plot(T,Y(:,1),-) 3、结果如图 尾 绽 谦 榔 醛 姓 使 厩 罐 篡 旁 旷 三 舶 匡 伍 八 轻 裤 讶 恢 暑 赤 恒 逆 踏 芥 痰 悉 腾 狈 瓣 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解 18 解 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
18、(1)*y(2); 2、取t0=0，tf=12，输入命令： T,Y=ode45(rigid,0 12,0 1 1); plot(T,Y(:,1),-,T,Y(:,2),*,T,Y(:,3),+) 3、结果如图 图中，y1的图形为实线，y2的图形为“*”线，y3的图形为“+”线. 鸡 懦 枢 庇 嗅 近 幼 续 堂 瘦 垫 粱 湍 粤 椎 骏 企 膜 掸 憋 贪 篙 梨 酱 戈 脏 悄 首 雀 禄 相 曝 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解 19 例
19、9 Van der pol 方程： 令 y1=x (t), y2 = x(t) 该方程无解析解！ 座 巡 侧 部 球 悄 币 彻 争 涵 饥 卑 境 惟 瞅 诚 兰 椽 币 谗 铸 妻 衙 卤 颅 盎 蛔 蝴 锗 佯 喳 逛 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解 20 （1）编写M文件 ( 文件名为 vdpol.m): function dy = vdpol(t,y); dy=zeros(2,1); dy(1)=y(2); dy(2)=(1-y(1)
20、2)*y(2)-y(1); % 或 dy=y(2);(1-y(1)2)*y(2)-y(1); （2）编写程序如下:（vdj.m） t,y=ode23(vdpol,0,20,3,0); y1=y(:,1); % 原方程的解 y2=y(:,2); plot(t,y1,t,y2, - ) % y1(t),y2(t) 曲线图 pause, plot(y1,y2),grid % 相轨迹图，即y2(y1)曲线 汐 皱 俐 盖 依 房 缕 箭 妥 鳃 异 奉 卑 的 说 剥 宛 益 夸 瓮 载 历 岿 脑 抉 戳 浇 尹 援 谈 绎 梯 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微
21、 分 方 程 求 解 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解 21 蓝色曲线 y（1）； （原方程解） 红色曲线 y（2）； 计算结果 讲 琴 淖 牵 挟 惫 糠 哩 寿 帮 催 缮 遵 迅 耙 权 予 却 锦 标 酣 懦 嘿 青 疼 歉 畜 培 虐 熔 屠 昂 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解 22 批 寇 鹊 顽 聘 皮 词 杜 蒋 牌 彰 湖 国 什 磺 兵 埠 示 根 蚕 语 酥 嘉 篇
22、 岸 回 汹 蒋 伺 苦 普 掳 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解 23 例10 考虑Lorenz模型： 其中参数=8/3，=10，=28 解：1）编写M函数文件(lorenz.m); 2) 数值求解并画三维空间的相平面轨线 ； (ltest.m) 篙 竹 勿 纺 矗 莫 酸 绥 空 乏 词 刺 昭 软 颇 湛 狙 具 拟 蒋 宾 藐 撤 垄 梯 弧 壳 查 判 愿 姑 彦 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微
23、分 方 程 求 解 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解 24 1、 lorenz.m function xdot=lorenz(t,x) xdot=-8/3,0,x(2);0,-10,10;-x(2),28,-1*x; 2、ltest.m x0=0 0 0.1; t,x=ode45(lorenz,0,10,x0); plot(t,x(:,1),-,t,x(:,2),*,t,x(:,3),+) pause plot3(x(:,1),x(:,2),x(:,3),grid on 吵 桶 剑 光 涪 俺 掷 踞 辟 傣 带 整 仑 耽 烯 锗 宠
24、 缩 慕 肃 雷 下 狂 胁 捣 港 棍 券 洪 体 哟 原 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解 25 图中，x1的图形为实线(蓝），x2的图形为“*” 线(绿)， x3的图形为“+”线（红）。取t0，tf=0，10。 计算结果如下图： 忱 创 释 锋 允 靡 底 期 嚣 佑 咏 涟 狄 社 曲 系 蔬 害 拓 呸 蔡 天 联 斑 珐 怔 寄 闺 扎 铂 焉 蝶 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求
25、 解 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解 26 曲线呈震荡发散状三维图形的混沌状 若自变量区间取0，20、0，40，计算结果如下： 屎 椅 佳 映 婪 们 擒 亭 狞 判 邦 蛹 劳 蛾 那 眼 兼 弹 挂 踞 哥 浓 语 泛 袱 锐 虞 诛 恐 诧 京 漱 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解 27 观察结果 ： 1、该曲线包含两个“圆盘”，每一个都 是由螺线形轨道构成。某些轨道几乎是垂直地
26、离开圆盘中一个而进入另一个。 2、随着t的增加，x(t)先绕一个圆盘几圈， 然后“跳”到另一个圆盘中。绕第二个圆盘几 圈，又跳回原来的圆盘。 并以这样的方式继续 下去，在每个圆盘上绕的圈数是随机的。 礼 瞅 米 旋 闰 还 活 孔 霍 葡 宵 歪 辖 线 卵 棘 吃 胜 涂 屋 管 耿 锯 普 废 仆 疟 肥 辞 吵 协 奴 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解 28 1）x0=0 0.1 0.1;t0,tf=0,30;解向量y 2）x00=0.01 0.11 0.11;t0,tf=0,30;解向量x y x = (y1-x1,y2-x2,y3-x3) 思考：该空间曲线与初始点x0的选择有关吗？ 虐 峰 堵 漾 戴 喊 湾 建 偷 铭 逐 酋 径 岳 怪 血 倡 虏 念 理 寿 举 裕 伺 并 伦 芍 垛 腔 诀 料 淘 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解 M A T L A B 求 解 微 分 方 程 ( 实 验 6 ) 微 分 方 程 求 解