实验二： 微分方程与差分方程模型Matlab 求解
一、实验目的
[1] 掌握解析、数值解法，并学会用图形观察解的形态和进行解的定性分析； [2] 熟悉MATLAB 软件关于微分方程求解的各种命令；
[3] 通过范例学习建立微分方程方面的数学模型以及求解全过程； [4] 熟悉离散 Logistic 模型的求解与混沌的产生过程。
二、实验原理
1. 微分方程模型与MATLAB 求解
解析解
用MATLAB 命令dsolve(‘eqn1’,’eqn2’, ...) 求常微分方程（组）的解析解。其中‘eqni'表示第i 个微分方程，Dny 表示y 的n 阶导数,默认的自变量为t 。 （1） 微分方程 例1 求解一阶微分方程 21y dx
dy
+= (1) 求通解 输入：
dsolve('Dy=1+y^2')
输出：
ans =
tan(t+C1)
（2）求特解 输入：
dsolve('Dy=1+y^2','y(0)=1','x')
指定初值为1，自变量为x 输出：
ans =
tan(x+1/4*pi)
例2 求解二阶微分方程 221
()0
4
(/2)2(/2)2/x y xy x y y y πππ
'''++-=='=-
原方程两边都除以2x ，得211
(1)04y y y x x
'''++-= 输入:
dsolve('D2y+(1/x)*Dy+(1-1/4/x^2)*y=0','y(pi/2)=2,Dy(pi/2)=-2/pi','x')
ans =
- (exp(x*i)*(pi/2)^(1/2)*i)/x^(1/2) +
(exp(x*i)*exp(-x*2*i)*(pi/2)^(3/2)*2*i)/(pi*x^(1/2))
试试能不用用simplify 函数化简 输入: simplify(ans)
ans =
2^(1/2)*pi^(1/2)/x^(1/2)*sin(x) （2）微分方程组
例3 求解 d f /d x =3f +4g ; d g /d x =-4f +3g 。
（1）通解:
[f,g]=dsolve('Df=3*f+4*g','Dg=-4*f+3*g')
f =
exp(3*t)*(C1*sin(4*t)+C2*cos(4*t)) g =
exp(3*t)*(C1*cos(4*t)-C2*sin(4*t))
特解：
[f,g]=dsolve('Df=3*f+4*g','Dg=-4*f+3*g','f(0)=0,g(0)=1')
f =
exp(3*t)*sin(4*t) g =
exp(3*t)*cos(4*t)
数值解
在微分方程(组)难以获得解析解的情况下，可以用Matlab 方便地求出数值解。格式为:
[t,y] = ode23('F',ts,y0,options)
注意：
? 微分方程的形式：y ' = F (t , y )，t 为自变量，y 为因变量（可以是多个，
如微分方程组）；
? [t, y]为输出矩阵，分别表示自变量和因变量的取值；
? F 代表一阶微分方程组的函数名（m 文件，必须返回一个列向量，每个
元素对应每个方程的右端）；
? ts 的取法有几种，（1）ts=[t0, tf] 表示自变量的取值范围，（2）
ts=[t0,t1,t2,…,tf],则输出在指定时刻t0,t1,t2,…,tf 处给出，（3）ts=t0:k:tf,则输出在区间[t0,tf]的等分点给出； ? y0为初值条件；
? options 用于设定误差限（缺省是设定相对误差是10^(-3)，绝对误差
是10^(-6)）； ode23是微分方程组数值解的低阶方法，ode45为中阶方法，与ode23类似。
例4 求解一个经典的范得波（Van Der pol ）微分方程：
0)0(',1)0(0')1(''2===+-+u u u u u u ，
解 形式转化：令)(');(21t u y t u y ==。则以上方程转化一阶微分方程组：
1221221
)1(;y y y y y y --='='。 编写M 文件如下，必须是M 文件表示微分方程组，并保存，一般地，M 文件的名字与函数名相同，保存位置可以为默认的work 子目录，也可以保存在自定义文件夹，这时注意要增加搜索路径（File\Set Path\Add Folder ）
function dot1=vdpol(t,y); dot1=[y(2); (1-y(1)^2)*y(2)-y(1)]; 在命令窗口写如下命令:
[t,y]=ode23('vdpol',[0,20],[1,0]); y1=y(:,1);y2=y(:,2);
plot(t,y1,t,y2,'--');title('Van Der Pol Solution ');
xlabel('Time,Second');ylabel('y(1)andy(2)')
注：Van der Pol方程描述具有一个非线性振动项的振动子的运动过程。最初,由于它在非线性电路上的应用而引起广泛兴趣。一般形式为
-
+u
+
u
u。
u
''2=
'
)1
(
图形解
无论是解析解还是数值解，都不如图形解直观明了。即使是在得到了解析解或数值解的情况下，作出解的图形，仍然是一件深受欢迎的事。这些都可以用Matlab方便地进行。
(1)图示解析解
如果微分方程（组）的解析解为：y=f (x)，则可以用Matlab函数fplot作出其图形：
fplot('fun',lims)
其中：fun给出函数表达式；lims=[xmin xmax ymin ymax]限定坐标轴的大小。例如
fplot('sin(1/x)', [0.01 0.1 -11])
(2)图示数值解
设想已经得到微分方程（组）的数值解(x,y)。可以用Matlab函数plot(x,y)直接作出图形。其中x和y为向量（或矩阵）。
2、Volterra模型（食饵捕食者模型）
意大利生物学家Ancona曾致力于鱼类种群相互制约关系的研究，他从第一次世界大战期间,地中海各港口捕获的几种鱼类捕获量百分比的资料中，发现鲨鱼的比例有明显增加（见下表）。
战争为什么使鲨鱼数量增加？是什么原因？
因为战争使捕鱼量下降，食用鱼增加，显然鲨鱼也随之增加。
但为何鲨鱼的比例大幅增加呢？生物学家Ancona无法解释这个现象，于是求助于著名的意大利数学家V.Volterra，希望建立一个食饵—捕食者系统的数学模型，定量地回答这个问题。
1、符号说明：
①x
1(t), x
2
(t)分别是食饵、捕食者（鲨鱼）在t时刻的数量；
②r1, r2是食饵、捕食者的固有增长率；
③λ1是捕食者掠取食饵的能力，
λ2是食饵对捕食者的供养能力；
2、基本假设：
① 捕食者的存在使食饵的增长率降低，假设降低的程度与捕食者数量成正比，
即)(21111
x r x dt
dx λ-=
② 食饵对捕食者的数量x 2起到增长的作用， 其程度与食饵数量x 1成正比，
即)(12222
x r x dt
dx λ+-=
综合以上①和②，得到如下模型：
模型一：不考虑人工捕获的情况
该模型反映了在没有人工捕获的自然环境中食饵与捕食者之间的制约关
系，没有考虑食饵和捕食者自身的阻滞作用，是Volterra 提出的最简单的模型。
给定一组具体数据，用matlab 软件求解。 食饵： r 1= 1, λ1= 0.1, x 10= 25; 捕食者(鲨鱼)：r 2=0.5, λ2=0.02, x 20= 2;
编制程序如下
1、建立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)); 2、在命令窗口执行如下程序：
[t,x]=ode45('shier',0:0.1:15,[25 2]); plot(t,x(:,1),'-',t,x(:,2),'*'),grid
??????
?+-=-=)
()(1222221111
x r x dt
dx x r x dt dx λλ2
)0(,25)0()
02.05.0()1.01(21122211
==??????
?+-=-=x x x x dt
dx x x dt dx
从图中可以看出它们的数量都呈现出周期性，而且鲨鱼数量的高峰期稍滞后于食饵数量的高峰期。
画出相轨迹图：
模型二 考虑人工捕获的情况
假设人工捕获能力系数为e ，相当于食饵的自然增长率由r 1 降为r 1-e ，捕食者的死亡率由r 2 增为 r 2+e ，因此模型一修改为：
设战前捕获能力系数e=0.3, 战争中降为e=0.1, 其它参数与模型(一)的参
?????++-=--=])([])[(1222221111
x e r x dt
dx
x e r x dt dx λλ
数相同。观察结果会如何变化？
1）当e=0.3时： 2）当e=0.1时：
分别求出两种情况下鲨鱼在鱼类中所占的比例。
即计算
画曲线：plot(t,p1(t),t,p2(t),'*') MATLAB 编程实现 建立两个M 文件
function dx=shier1(t,x) dx=zeros(2,1);
dx(1)=x(1)*(0.7-0.1*x(2)); dx(2)=x(2)*(-0.8+0.02*x(1)); function dy=shier2(t,y) dy=zeros(2,1);
dy(1)=y(1)*(0.9-0.1*y(2)); dy(2)=y(2)*(-0.6+0.02*y(1));
运行以下程序：
[t1,x]=ode45('shier1',[0 15],[25 2]); [t2,y]=ode45('shier2',[0 15],[25 2]); x1=x(:,1);x2=x(:,2); p1=x2./(x1+x2);
y1=y(:,1);y2=y(:,2); p2=y2./(y1+y2);
1
122
2112(0.70.1)(0.80.02)(0)25,(0)2dx x x dt dx x x dt
x x ?=-???=-+??==???
????
?????==+-=-=2)0(,25)0()02.06.0()1.09.0(21122
211
y y y y dt
dy y y dt dy )
()()
()(;
)
()()
()(21222121t y t y t y t p t x t x t x t p +=
+=
plot(t1,p1,'-',t2,p2,'*')
051015
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
图中‘*’曲线为战争中鲨鱼所占比例。 结论：战争中鲨鱼的比例比战前高。
3、 Logistic 映射
logistic 映射---通向混沌的道路
混沌系统，由于其行为的复杂性，往往认为其动态特性（运动方程）也一定非常复杂，事实并非如此，一个参量很少、动态特性非常简单的系统有时也能够产生混沌现象，以一维虫口模型为例，假设某一区域内的现有虫口数为y n ，昆虫的繁殖率为r ，且第n 代昆虫不能存活于第n+1代，既无世代交叠，则第n+1代虫口数为1n n y ry +=，r ＞1时，虫口会无限制地增长；r ＜1时，虫口最终会趋于消亡，因此需要对模型进行修正。由于环境的制约和食
物有限，因争夺生存空间发生相互咬斗事件的最大次数为
1
(1)2
n n y y -，即制约虫口数的因素与2n y 成正比，设咬斗事件的战死率为β则对虫口的修正项为2
n y β- ，则有：
2
1n n n y ry y β+=-.
令n n x y r
β
=
，则
1(1)n n n x rx x +=- （1）
取最大虫口数为1，且虫口数不能为负，则 ；当 =0.5时，方程有极大
值，11
4
n x r +=
而 又必须小于1，因而r ＜4，则参量r 的取值范围为1到4，这就得到
一个抽象的标准虫口方程（1）。
记映射f 为
()(1)f x rx x =- （2）
方程（1）可写为
1()n n x f x += （3）
这一迭代关系通常称为logistic 映射。从[0,1]内点x 0出发，由Logistic 映射的
迭代形成
x n = f n
(x 0), n = 0,1,2,…
序列{x n }称为x 0的轨道。
一个看似简单的系统，随着参量的不同会表现出截然不同的行为，当r 的取值范围
在1~3时，方程（1）有定态解
1
1x r
*=-
即方程通过多次迭代后趋于一个稳定的不动点，此时系统是稳定的。
1
1x r
*=-为方程()f x x =的解，称为周期2点。
当r 在3~3.448范围内取值时，经过多次迭代，方程（1）出现周期2点1x *
和2x *
，即
1x *和2x *是方程2
()f x x =的解，满足
122(1)x rx x ***
=-
211(1)x rx x ***=-
3r =是使解有意义的r 最小值。
随着r 的增大，r=3.449；3.544；3.564…依次出现周期4、周期8、周期16…的振
荡解，r 的极限值约为3.569。这种行为称为倍周期分岔，直到r ＞3.5699时，系统进入了混沌状态，如下图所示，此时系统的状态不再具有规律性，而是发生随机的波动，使图d 的右侧的大部分区域被涂黑了，仔细观察发现，混沌区域并非一片涂斑，而是有粗粗细细的
白色“竖线”，称为周期窗口，随着参量r 的增大（如 1r =+系统出现周期三的稳定状态，
Logistic映射分岔图
接着倍周期分岔以更快的速度进行，再次进入混沌状态。如果将周期窗口放大，发现其结构与分岔图的整体结构具有相似性，而且是一种无限嵌套的自相似结构。
可以看出，通过改变系统参量，使系统进入混沌的第一种模式是倍周期分岔，即由不动点→周期二→周期四→…无限倍周期→进入混沌状态。当然通向混沌的道路不只于此，第二种通向的道路是：从平衡态到周期运动，再到拟周期运动，直到进入混沌状态。第三种通向混沌的方式是阵发（或间歇）道路，即系统在近似周期运动的过程中，通过改变参量，系统会出现阵发性混沌过程，随着参量的调整，阵发性混沌越来越频繁，近似的周期运动越来越少，最后进入混沌。
Matlab程序
下面程序绘制r=2，x0=0.3的轨道
clear all;clf;
x=0.3;
r=2;
n=input('Please input a number: ');
for i=1:n
x1=r*x*(1-x);
x=x1;
plot(i,x1,'.');
hold on;
end
下面程序绘制logistic映射r=3,5的轨道，观察
是否有周期点
clear all;clf;
x=0.3;r=3.5;
for i=1:100
x1=r*x*(1-x);
x=x1;
plot(i,x1,'.'); hold on ; end
下面程序绘制logistic 映射r=4的轨道，观察其混沌
clear all ;clf; x=0.007; for i=1:100
x1=4*x*(1-x); x=x1;
plot(i,x1,'.'); hold on ; end
下面程序绘制在混沌状态下不同初值x 0=0.100001和x 0=0.1的轨道的差（对初值的敏感性）
clear all ;clf;
x1=0.100001;x11=0.1; for i=1:1000
x2=4*x1*(1-x1); x1=x2;
x22=4*x11*(1-x11); x11=x22;
plot(i,x11-x1); hold on ; end
下面程序绘制logistic 映射的分岔图
clear all ;clf; x=0.2;
for r=1:0.001:4 for i=1:18
x1=r*x*(1-x); x=x1; if i>9
plot(r,x); hold on ; end end end
1 1.5
2 2.5
3 3.54
蛛网迭代
clear all;clf;
x=0.007;y=0;
for i=1:200
XY=[x,y];
y=4*x*(1-x);
XY1=[x,y];
line(XY,XY1);hold on;
x=y;
XY2=[x,x];
line(XY1,XY2);hold on;
end
三、实验内容
1．求微分方程的解析解, 并画出图形， =+2,(0)1,01y y x y x '=<< 2．求微分方程的数值解, 并画出图形， cos 0,(0)1,(0)0y y x y y '''+===
3．两种相似的群体之间为了争夺有限的同一种食物来源和生活空间而进行生存竞争时，往往是竞争力较弱的种群灭亡，而竞争力较强的种群达到环境容许的最大数量。
假设有甲、乙两个生物种群，当它们各自生存于一个自然环境中，均服从 Logistic 规律。
（1）)(),(21t x t x 是两个种群的数量； （2）21,r r 是它们的固有增长率； （3）21,n n 是它们的最大容量；
（4）)(12m m 为种群乙(甲)占据甲(乙)的位置的数量，并且 .1122;x m x m βα==
计算)(),(21t x t x , 画出图形及相轨迹图。解释其解变化过程。
2）121,r r ==设12100,n n ==102010x x ==,α=1.5，β=0.7，计算)(),(21t x t x ,
画出图形及相轨迹图。解释其解变化过程。
4. 绘制Logistic 映射的轨道图,分岔图和蛛网迭代图。
四、实验心得
利用MATLAB求解常微分方程数值解
目录 1. 内容简介 (1) 2. Euler Method(欧拉法)求解 (1) 2.1. 显式Euler法和隐式Euler法 (2) 2.2. 梯形公式和改进Euler法 (3) 2.3. Euler法实用性 (4) 3. Runge-Kutta Method(龙格库塔法)求解 (5) 3.1. Runge-Kutta基本原理 (5) 3.2. MATLAB中使用Runge-Kutta法的函数 (7) 4. 使用MATLAB求解常微分方程 (7) 4.1. 使用ode45函数求解非刚性常微分方程 (8) 4.2. 刚性常微分方程 (9) 5. 总结 (9) 参考文献 (11) 附录 (12) 1. 显式Euler法数值求解 (12) 2. 改进Euler法数值求解 (12) 3. 四阶四级Runge-Kutta法数值求解 (13) 4.使用ode45求解 (14)
1.内容简介 把《高等工程数学》看了一遍，增加对数学内容的了解，对其中数值解法比较感兴趣，这大概是因为在其它各方面的学习和研究中经常会遇到数值解法的问题。理解模型然后列出微分方程，却对着方程无从下手，无法得出精确结果实在是让人难受的一件事情。 实际问题中更多遇到的是利用数值法求解偏微分方程问题，但考虑到先从常微分方程下手更为简单有效率，所以本文只研究常微分方程的数值解法。把一个工程实际问题弄出精确结果远比弄清楚各种细枝末节更有意思，因此文章中不追求非常严格地证明，而是偏向如何利用工具实际求解出常微分方程的数值解，力求将课程上所学的知识真正地运用到实际方程的求解中去，在以后遇到微分方程的时候能够熟练运用MATLAB得到能够在工程上运用的结果。 文中求解过程中用到MATLAB进行数值求解，主要目的是弄清楚各个函数本质上是如何对常微分方程进行求解的，对各种方法进行MATLAB编程求解，并将求得的数值解与精确解对比，其中源程序在附录中。最后考察MATLAB中各个函数的适用范围，当遇到实际工程问题时能够正确地得到问题的数值解。 2.Euler Method(欧拉法)求解 Euler法求解常微分方程主要包括3种形式，即显式Euler法、隐式Euler法、梯形公式法，本节内容分别介绍这3种方法的具体内容，并在最后对3种方法精度进行对比，讨论Euler法的实用性。 本节考虑实际初值问题 使用解析法，对方程两边同乘以得到下式
用Matlab 软件求解微分方程 1．解析解 （1）一阶微分方程 求21y dx dy +=的通解：dsolve('Dy=1+y^2','x') 求y x dx dy -+=21的通解：dsolve('Dy=1+x^2-y','x') 求?????=+=1 )0(12y y dx dy 的特解：dsolve('Dy=1+y^2',’y(0)=1’,'x') （2）高阶微分方程 求解???-='==-+'+''. 2)2(,2)2(,0)(222πππy y y n x y x y x 其中，21=n ，命令为： dsolve('x^2*D2y+x*Dy+(x^2-0.5^2)*y=0','y(pi/2)=2,Dy(pi/2)=-2/pi','x') 求042=-+'-'''x y y y 的通解，命令为： dsolve('D3y-2*Dy+y-4*x=0','x') 输出为： ans=8+4*x+C1*exp(x)+C2*exp(-1/2*(5^(1/2)+1)*x)+C3*exp(1/2*(5^(1/2)-1)*x) （3）一阶微分方程组 求???+-='+='). (3)(4)(),(4)(3)(x g x f x g x g x f x f 的通解：[f,g]=dsolve('Df=3*f+4*g','Dg=-4*f+3*g','x') 输出为： f =exp(3*x)*(cos(4*x)*C1+sin(4*x)*C2) g =-exp(3*x)*(sin(4*x)*C1-cos(4*x)*C2) 若再加上初始条件1)0(,0)0(==g f ，则求特解： [f,g]=dsolve('Df=3*f+4*g','Dg=-4*f+3*g','f(0)=0,g(0)=1','x') 输出为： f =exp(3*x)*sin(4*x) g =exp(3*x)*cos(4*x) 2．数值解 （1）一阶微分方程
实验四 求微分方程的解 一、问题背景与实验目的 实际应用问题通过数学建模所归纳而得到的方程，绝大多数都是微分方程，真正能得到代数方程的机会很少．另一方面，能够求解的微分方程也是十分有限的，特别是高阶方程和偏微分方程（组）．这就要求我们必须研究微分方程（组）的解法，既要研究微分方程（组）的解析解法（精确解），更要研究微分方程（组）的数值解法（近似解）． 对微分方程（组）的解析解法(精确解)，Matlab 有专门的函数可以用，本实验将作一定的介绍． 本实验将主要研究微分方程(组)的数值解法（近似解），重点介绍 Euler 折线法． 二、相关函数（命令）及简介 1．dsolve ('equ1','equ2',…)：Matlab 求微分方程的解析解．equ1、equ2、…为方程（或条件）．写方程（或条件）时用 Dy 表示y 关于自变量的一阶导数，用用 D2y 表示 y 关于自变量的二阶导数，依此类推． 2．simplify(s )：对表达式 s 使用 maple 的化简规则进行化简． 例如： syms x simplify(sin(x)^2 + cos(x)^2) ans=1 3．[r,how]=simple(s)：由于 Matlab 提供了多种化简规则，simple 命令就是对表达式 s 用各种规则进行化简，然后用 r 返回最简形式，how 返回形成这种形式所用的规则． 例如： syms x [r,how]=simple(cos(x)^2-sin(x)^2) r = cos(2*x) how = combine 4．[T,Y] = solver(odefun,tspan,y 0) 求微分方程的数值解． 说明： (1) 其中的 solver 为命令 ode45、ode23、ode113、ode15s 、ode23s 、ode23t 、ode23tb 之一． (2) odefun 是显式常微分方程：?????==0 0)() ,(y t y y t f dt dy (3) 在积分区间 tspan =],[0f t t 上，从0t 到f t ，用初始条件0y 求解．
第四讲Matlab求解微分方程(组) 理论介绍:Matlab求解微分方程(组)命令 求解实例:Matlab求解微分方程(组)实例 实际应用问题通过数学建模所归纳得到得方程,绝大多数都就是微分方程,真正能得到代数方程得机会很少、另一方面,能够求解得微分方程也就是十分有限得,特别就是高阶方程与偏微分方程(组)、这就要求我们必须研究微分方程(组)得解法:解析解法与数值解法、 一.相关函数、命令及简介 1、在Matlab中,用大写字母D表示导数,Dy表示y关于自变量得一阶导数,D2y 表示y关于自变量得二阶导数,依此类推、函数dsolve用来解决常微分方程(组)得求解问题,调用格式为: X=dsolve(‘eqn1’,’eqn2’,…) 函数dsolve用来解符号常微分方程、方程组,如果没有初始条件,则求出通解,如果有初始条件,则求出特解、 注意,系统缺省得自变量为t 2、函数dsolve求解得就是常微分方程得精确解法,也称为常微分方程得符号解、但就是,有大量得常微分方程虽然从理论上讲,其解就是存在得,但我们却无法求出其解析解,此时,我们需要寻求方程得数值解,在求常微分方程数值解方 面,MATLAB具有丰富得函数,我们将其统称为solver,其一般格式为: [T,Y]=solver(odefun,tspan,y0) 说明:(1)solver为命令ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb、ode15i之一、 (2)odefun就是显示微分方程在积分区间tspan上从到用初始条件求解、 (3)如果要获得微分方程问题在其她指定时间点上得解,则令tspan(要求就是单调得)、 (4)因为没有一种算法可以有效得解决所有得ODE问题,为此,Matlab提供了多种求解器solver,对于不同得ODE问题,采用不同得solver、 表1 Matlab中文本文件读写函数
《高等数学》实验报告 实验人员:系(班): 学号: 姓名: 实验地点:电教楼五号机房 实验名称:Matlab 高等数学实验 实验时间:2014-6-3 16:30--18:30 实验名称:用Matlab 软件求常微分方程的解(或通解) 实验目的:熟练掌握Matlab 软件求常微分方程的解(或通解) 实验内容:(给出实验程序与运行结果) 1、求微分方程的特解. 1、?? ?? ?===+-10)0(,6)0(034'2 2y y y dx dy dx y d 程序:>> dsolve('D2y-4*Dy+3*y','y(0)=6,Dy(0)=10','x') ans = 4*exp(x)+2*exp(3*x) 吕梁学院《高等数学》实验报告 情况试高中
2、?? ???===++0)0(,2)0(044'2 2y y y dx dy dx y d 程序:>>dsolve('4*D2y+4*Dy+y','y(0)=2,Dy(0)=0','x') ans = 2*exp(-1/2*x)+exp(-1/2*x)*x 3、?? ???===++15)0(',0)0(029422y y y dx dy dx y d 程序:>>dsolve('D2y+4*Dy+29*y=0','y(0)=9,Dy(0)=15','x') ans = 33/5*exp(-2*x)*sin(5*x)+9*exp(-2*x)*cos(5*x) 4、?? ???===+-3)0(',0)0(013422y y y dx dy dx y d 程序:>>dsolve('D2y-4*dy+13*y=0','y(0)=0','Dy(0)=3','x') ans = 3/13*sin(13^(1/2)*x)*13^(1/2)-4/13*cos(13^(1/2)*x)*dy+4/13*dy 5、?? ???-===--5)0(',0)0(04322y y y dx dy dx y d 程序:>>dsolve('D2y-3*Dy-4*y','y(0)=0,Dy(0)=-5','x') ans = exp(-x)-exp(4*x)
第四讲 Matlab 求解微分方程（组） 理论介绍：Matlab 求解微分方程（组）命令 求解实例：Matlab 求解微分方程（组）实例 实际应用问题通过数学建模所归纳得到的方程，绝大多数都是微分方程，真正能得到代数方程的机会很少.另一方面，能够求解的微分方程也是十分有限的，特别是高阶方程和偏微分方程（组）.这就要求我们必须研究微分方程（组）的解法：解析解法和数值解法. 一．相关函数、命令及简介 1.在Matlab 中，用大写字母D 表示导数，Dy 表示y 关于自变量的一阶导数，D2y 表示y 关于自变量的二阶导数，依此类推.函数dsolve 用来解决常微分方程（组）的求解问题，调用格式为： X=dsolve(‘eqn1’,’eqn2’,…) 函数dsolve 用来解符号常微分方程、方程组，如果没有初始条件，则求出通解，如果有初始条件，则求出特解. 注意，系统缺省的自变量为t 2.函数dsolve 求解的是常微分方程的精确解法，也称为常微分方程的符号解.但是，有大量的常微分方程虽然从理论上讲，其解是存在的，但我们却无法求出其解析解，此时，我们需要寻求方程的数值解，在求常微分方程数值解方面，MATLAB 具有丰富的函数，我们将其统称为solver ，其一般格式为： [T,Y]=solver(odefun,tspan,y0) 说明：(1)solver 为命令ode45、ode23、ode113、ode15s 、ode23s 、ode23t 、ode23tb 、ode15i 之一. (2)odefun 是显示微分方程'(,)y f t y =在积分区间tspan 0[,]f t t =上从0t 到f t 用初始条件0y 求解. (3)如果要获得微分方程问题在其他指定时间点012,,, ,f t t t t 上的解，则令 tspan 012[,,,]f t t t t =(要求是单调的). (4)因为没有一种算法可以有效的解决所有的ODE 问题，为此，Matlab 提供
用Matlab件求常微分方程解(或通解)
————————————————————————————————作者：————————————————————————————————日期：
《高等数学》实验报告 实验人员:系(班): 学号: 姓名: 实验地点:电教楼五号机房 实验名称:Matlab 高等数学实验 实验时间:2014-6-3 16:30--18:30 实验名称:用Matlab 软件求常微分方程的解(或通解) 实验目的:熟练掌握Matlab 软件求常微分方程的解(或通解) 实验内容:(给出实验程序与运行结果) 一、求微分方程的特解. 1、?? ???===+-10)0(,6)0(034'22y y y dx dy dx y d 程序:>> dsolve('D2y-4*Dy+3*y','y(0)=6,Dy(0)=10','x') ans = 4*exp(x)+2*exp(3*x) 吕梁学院《高等数学》实验报告
2、?? ???===++0)0(,2)0(044'22y y y dx dy dx y d 程序:>>dsolve('4*D2y+4*Dy+y','y(0)=2,Dy(0)=0','x') ans = 2*exp(-1/2*x)+exp(-1/2*x)*x 3、?? ???===++15)0(',0)0(029422y y y dx dy dx y d 程序:>>dsolve('D2y+4*Dy+29*y=0','y(0)=9,Dy(0)=15','x') ans = 33/5*exp(-2*x)*sin(5*x)+9*exp(-2*x)*cos(5*x) 4、?? ???===+-3)0(',0)0(013422y y y dx dy dx y d 程序:>>dsolve('D2y-4*dy+13*y=0','y(0)=0','Dy(0)=3','x') ans = 3/13*sin(13^(1/2)*x)*13^(1/2)-4/13*cos(13^(1/2)*x)*dy+4/13*dy 5、?? ???-===--5)0(',0)0(04322y y y dx dy dx y d 程序:>>dsolve('D2y-3*Dy-4*y','y(0)=0,Dy(0)=-5','x') ans = exp(-x)-exp(4*x)
数学应用软件作业6-用Matlab 求解微分方程(组)的解析解和数值解
注：上机作业文件夹以自己的班级姓名学号命名，文件夹包括如下上机报告和Matlab程序。 上机报告模板如下： 佛山科学技术学院 上机报告 课程名称数学应用软件 上机项目用Matlab求解微分方程(组)的解析解和数值解 专业班级姓名学号 一. 上机目的 1.了解求微分方程（组）的解的知识。 2.学习Matlab中求微分方程的各种解的函数，如 dsolve命令、ode45函数等等，其中注意把方程化为新的方程的形式。 3.掌握用matlab编写程序解决求解微分方程的问 题。 二. 上机内容 1、求高阶线性齐次方程：y’’’-y’’-3y’+2y=0。 2、求常微分方程组
2 210cos,2 24,0 t t t dx dy x t x dt dt dx dy y e y dt dt = - = ? +-== ?? ? ?++== ?? 3、求解 分别用函数ode45和ode15s计算求解，分别画出图形，图形分别标注标题。 4、求解微分方程 ,1 )0( ,1 '= + + - =y t y y 先求解析解，在[0,1]上作图； 再用ode45求数值解(作图的图形用“o”表示)，在同一副图中作图进行比较，用不同的颜色表示。 三. 上机方法与步骤 给出相应的问题分析及求解方法，并写出Matlab 程序，并有上机程序显示截图。 题1：直接用命令dsolve求解出微分方程的通解。 Matlab程序：
dsolve('D3y-D2y-3*Dy+2*y','x') 题2：将微分方程组改写为 5cos2exp(2) 5cos2exp(2) (0)2,(0)0 dx t t x y xt dy t t x y dt x y ? =+--- ? ? ? =-+-+- ? ? == ? ? ? ， 再用命令dsolve求解微分方程的通解。 Matlab程序： 建立timu2.m如下： [x,y]=dsolve('Dx=5*cos(t)+2*exp(-2*t)-x-y','Dy=-5*cos(t)+2*exp(-2*t)+x-y ','x(0)=2,y(0)=0','t') x=simple(x) y=simple(y)
数学模型之微分方程及其MATLAB求解 ---卫星轨迹等经典例题求解分析1. 考虑初值问题画图 y'''?3y ''?y 'y = 0 y(0) = 0 y '(0) =1 y ' '(0) = ?1 2、 3、 【实验步骤与程序】 1. M -文件建立m函数文件
function y=f(t,x) y=[x(2);x(3);9*x(3)^2+x(1)*x(2)]; 求解微分方程，命令如下： x0=[0;1;-1]; [t,y]=ode45(@mm,[0,2.5],x0); plot(y(:,1),y(:,2)); figure(2); plot3(y(:,1),y(:,2),y(:,3))
2、M -文件建立m函数文件 function dx=appollo(t,x) mu=1/82.45; mustar=1-mu; r1=sqrt((x(1)+mu)^2+x(3)^2); r2=sqrt((x(1)-mustar)^2+x(3)^2); dx=[x(2) 2*x(4)+x(1)-mustar*(x(1)+mu)/r1^3-mu*(x(1)-mustar)/r2^3 x(4) -2*x(2)+x(3)-mustar*x(3)/r1^3-mu*x(3)/r2^3];
求解微分方程，命令如下： x0=[1.2;0;0;-1.04935751]; options=odeset('reltol',1e-8); [t,y]=ode45(@appollo,[0,20],x0,options); plot(y(:,1),y(:,3)) title('Appollo卫星运动轨迹') xlabel('x') ylabel('y')
MATLAB简介 MATLAB是矩阵实验室（Matrix Laboratory）的简称，是美国MathWorks公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境。 一、基本功能 MATLAB是由美国mathworks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式，代表了当今国际科学计算软件的先进水平。 MATLAB和Mathematica、Maple并称为三大数学软件。它在数学类科技应用软件中在数值计算方面首屈一指。MATLAB可以进行矩阵运算、绘制函数和数据、实现算法、创建用户界面、连接其他编程语言的程序等，主要应用于工程计算、控制设计、信号处理与通讯、图像处理、信号检测、金融建模设计与分析等领域。 MATLAB的基本数据单位是矩阵，它的指令表达式与数学、工程中常用的形式十分相似，故用MATLAB来解算问题要比用C，FORTRAN等语言完成相同的事情简捷得多，并且MATLAB也吸收了像Maple等软件的优点，使MATLAB成为一个强大的数学软件。在新的版本中也加入了对C，FORTRAN，C++，JAVA的支持。可以直接调用,用户也可以将自己编写的实用程序导入到MATLAB函数库中方便自己以后调用，此外许多的MATLAB爱好者都编写了一些经典的程序，用户可以直接进行下载就可以用。 二、特点 1) 高效的数值计算及符号计算功能，能使用户从繁杂的数学运算分析中解脱出来； 2) 具有完备的图形处理功能,实现计算结果和编程的可视化; 3) 友好的用户界面及接近数学表达式的自然化语言，使学者易于学习和掌握； 4) 功能丰富的应用工具箱(如信号处理工具箱、通信工具箱等) ,为用户提供了大量方便实用的处理工具。 三、优势 1.友好的工作平台编程环境 MATLAB由一系列工具组成。这些工具方便用户使用MATLAB的函数和文件，其中许多工具采用的是图形用户界面。包括MATLAB桌面和命令窗口、历史命令窗口、编辑器和调试器、路径搜索和用于用户浏览帮助、工作空间、文件的浏览器。随着MATLAB的商业化以及软件本身的不断升级，MATLAB的用户界面也越来越精致，更加接近Windows的标准界面，人机交互性更强，操作更简单。 2.强大的科学计算机数据处理能力 MATLAB是一个包含大量计算算法的集合。其拥有600多个工程中要用到的数学运算函数，可以方便的实现用户所需的各种计算功能，可以用它来代替底层编程语言，如C和C++ 。在计算要求相同的情况下，使用MATLAB的编程工作量会大大减少。
利用matlab编写S函数求解微分方程自动化专业综合设计报告 自动化专业综合设计报告
函数求解微S编写设计题目：利用 matlab 分方程 自动化系统仿真实验室所在 实验室： 郭卫平 指导教师： 律迪迪学生姓名 200990519114 班级文自0921 学号 成绩评定： 自动化专业综合设计报告
一、设计目的 了解使用simulink的扩展工具——S-函数，s函数可以利用matlab的丰富资源，而不仅仅局限于simulink提供的模块，而用c或c++等语言写的s函数还可以实现对硬件端口的操作，还可以操作windows API 等的，它的魅力在于完美结合了simulink 框图简洁明快的特点和编程灵活方便的优点，提供了增强和扩展sinulink能力的强大机制，同时也是使用RTW实现实时仿真的关键。 二、设计要求 求解解微分方程 y'=y-2x/y 自动化专业综合设计报告 y(0)=1 要求利用matlab编写S函数求解 三、设计内容（可加附页） 【步骤1】获取状态空间表达式。
在matlab中输入 dsolve(‘Dy=y-2*x/y','y(0)=1'， 'x') 得到 y=(2*x+1).^(1/2); 【步骤2】建立s函数的m文件。 利用21·用S函数模板文件。 以下是修改之后的模板文件sfuntmpl.m 的内容。 function [sys,x0,str,ts] = sfuntmpl(t,x,u,flag) %SFUNTMPL S-function M-file General template define you can With % M-file S-functions, you own ordinary differential system equations (ODEs), discrete % equations, and/or just about any type of algorithm to be used within a %
实验六 用matlab 求解常微分方程 1．微分方程的概念 未知的函数以及它的某些阶的导数连同自变量都由一已知方程联系在一起的方程称为微分方程。如果未知函数是一元函数，称为常微分方程。常微分方程的一般形式为 0),,",',,()(=n y y y y t F 如果未知函数是多元函数，成为偏微分方程。联系一些未知函数的一组微分方程组称为微分方程组。微分方程中出现的未知函数的导数的最高阶解数称为微分方程的阶。若方程中未知函数及其各阶导数都是一次的，称为线性常微分方程，一般表示为 )()(')()(1)1(1)(t b y t a y t a y t a y n n n n =++++-- 若上式中的系数n i t a i ,,2,1),( =均与t 无关，称之为常系数。 2．常微分方程的解析解 有些微分方程可直接通过积分求解.例如,一解常系数常微分方程1+=y dt dy 可化为 dt y dy =+1,两边积分可得通解为 1-=t ce y .其中c 为任意常数.有些常微分方程可用一些技巧,如分离变量法,积分因子法,常数变异法,降阶法等可化为可积分的方程而求得解析解. 线性常微分方程的解满足叠加原理,从而他们的求解可归结为求一个特解和相应齐次微分方程的通解.一阶变系数线性微分方程总可用这一思路求得显式解。高阶线性常系数微分方程可用特征根法求得相应齐次微分方程的基本解，再用常数变异法求特解。 一阶常微分方程与高阶微分方程可以互化，已给一个n 阶方程 ),,",',()1()(-=n n y y y t f y 设)1(21,,',-===n n y y y y y y ，可将上式化为一阶方程组 ?????????====-),,,,(''''2113221n n n n y y y t f y y y y y y y 反过来，在许多情况下，一阶微分方程组也可化为高阶方程。所以一阶微分方程组与高阶常微分方程的理论与方法在许多方面是相通的，一阶常系数线性微分方程组也可用特征根法求解。 3．微分方程的数值解法 除常系数线性微分方程可用特征根法求解，少数特殊方程可用初等积分法求解外，大部分微分方程无限世界，应用中主要依靠数值解法。考虑一阶常微分方程初值问题 ???=<<=000)()),(,()('y t y t t t t y t f t y f
Euler法解常微分方程 Euler法解常微分方程算法： Step 1 分别取积分上限、积分下限、步长 Step 2计算判断是否成立，成立转到Step 3，否则继续进行Step 4 Step 3 计算 Step 4 Euler法解常微分方程算程序： function euler2(fun,y0,A,h) %fun--y' %y0---初值 %A----x取值范围 %a----x左区间端点值 %b----x右区间端点值 %h----给定步长 x=min(A); b=max(A); y=y0; while x Step 3 （1）做显性Euler预测 （2）将带入 Step 4计算判断是否成立，成立返回Step 3，否则继续进行Step 5 Step 5 改进Euler法解常微分方程算程序： function gaijineuler2(fun,y0,A,h) %fun--y' %y0---初值 %A----x取值范围 %a----x左区间端点值 %b----x右区间端点值 %h----给定步长 a=min(A); b=max(A); x=a:h:b; y(1)=y0; for i=1:length(x)-1 w1=feval(fun,x(i),y(i)); y(i+1)=y(i)+h*w1; w2=feval(fun,x(i+1),y(i+1)); y(i+1)=y(i)+h*(w1+w2)/2; end x=x' y=y' 例：用改进Euler法计算下列初值问题（取步长h=0.25） 输入：fun=inline('-x*y^2') gaijineuler2(fun,2,[0 5],0.25) 得到： x = 0.2500 0.5000 0.7500 1.0000 1.2500 1.5000 1.7500 2.0000 2.2500 2.5000 2.7500 第四讲 Matlab 求解微分方程（组） 理论介绍：Matlab 求解微分方程（组）命令 求解实例：Matlab 求解微分方程（组）实例 实际应用问题通过数学建模所归纳得到的方程，绝大多数都是微分方程，真正能得到代数方程的机会很少.另一方面，能够求解的微分方程也是十分有限的，特别是高阶方程和偏微分方程（组）.这就要求我们必须研究微分方程（组）的解法：解析解法和数值解法. 一．相关函数、命令及简介 1.在Matlab 中，用大写字母D 表示导数，Dy 表示y 关于自变量的一阶导数，D2y 表示y 关于自变量的二阶导数，依此类推.函数dsolve 用来解决常微分方程（组）的求解问题，调用格式为： X=dsolve(‘eqn1’,’eqn2’,…) 函数dsolve 用来解符号常微分方程、方程组，如果没有初始条件，则求出通解，如果有初始条件，则求出特解. 注意，系统缺省的自变量为t 2.函数dsolve 求解的是常微分方程的精确解法，也称为常微分方程的符号解.但是，有大量的常微分方程虽然从理论上讲，其解是存在的，但我们却无法求出其解析解，此时，我们需要寻求方程的数值解，在求常微分方程数值解方面，MATLAB 具有丰富的函数，我们将其统称为solver ，其一般格式为： [T,Y]=solver(odefun,tspan,y0) 说明：(1)solver 为命令ode45、ode23、ode113、ode15s 、ode23s 、ode23t 、ode23tb 、ode15i 之一. (2)odefun 是显示微分方程'(,)y f t y =在积分区间tspan 0[,]f t t =上从0t 到f t 用初始条件0y 求解. (3)如果要获得微分方程问题在其他指定时间点012,,,,f t t t t 上的解，则令 tspan 012[,,, ]f t t t t =(要单调的). (4)因为没有一种算法可以有效的解决所有的ODE 问题，为此，Matlab 提供了多种求解器solver ，对于不同的ODE 问题，采用不同的solver. 常微分方程： 1 ODE解算器简介(ode**) 2 微分方程转换 3 刚性/非刚性问题(Stiff/Nonstiff) 4 隐式微分方程(IDE) 5 微分代数方程(DAE) 6 延迟微分方程(DDE) 7 边值问题(BVP) 偏微分方程(PDEs)Matlab解法 偏微分方程： 1 一般偏微分方程组(PDEs)的命令行求解 2 特殊偏微分方程(PDEs)的PDEtool求解 3 陆君安《偏微分方程的MATLAB解法 先来认识下常微分方程(ODE)初值问题解算器(solver) [T,Y,TE,YE,IE] = odesolver(odefun,tspan,y0,options) sxint = deval(sol,xint) Matlab中提供了以下解算器： 输入参数： odefun：微分方程的Matlab语言描述函数，必须是函数句柄或者字符串，必须写成Matlab 规范格式(也就是一阶显示微分方程组)，这个具体在后面讲解 tspan=[t0 tf]或者[t0,t1,…tf]：微分变量的范围，两者都是根据t0和tf的值自动选择步长，只是前者返回所有计算点的微分值，而后者只返回指定的点的微分值，一定要注意对于后者tspan必须严格单调，还有就是两者数据存储时使用的内存不同(明显前者多)，其它没有任何本质的区别 y0=[y(0),y’(0),y’’(0)…]：微分方程初值，依次输入所有状态变量的初值，什么是状态变量在后面有介绍 options：微分优化参数，是一个结构体，使用odeset可以设置具体参数，详细内容查看帮助 输出参数： T：时间列向量，也就是ode**计算微分方程的值的点 Y：二维数组，第i列表示第i个状态变量的值，行数与T一致 在求解ODE时，我们还会用到deval()函数，deval的作用就是通过结构体solution计算t 对应x值，和polyval之类的很相似！ 参数格式如下： sol：就是上次调用ode**函数得道的结构体解 xint：需要计算的点，可以是标量或者向量，但是必须在tspan范围内 该函数的好处就是如果我想知道t=t0时的y值，不需要重新使用ode计算，而直接使用上次计算的得道solution就可以 [教程] 微分方程转换为一阶显示微分方程组方法 好，上面我们把Matlab中的常微分方程(ODE)的解算器讲解的差不多了，下面我们就具体开始介绍如何使用上面的知识吧！ 现实总是残酷的，要得到就必须先付出，不可能所有的ODE一拿来就可以直接使用，因此，在使用ODE解算器之前，我们需要做的第一步，也是最重要的一步，借助状态变量将微分 Euler 法解常微分方程 Euler 法解常微分方程算法： Step 1 分别取积分上限、积分下限、步长 Step 2计算h n n +=判断b n ≤是否成立，成立转到Step 3，否则继续进行Step 4 Step 3 计算),(1n n n n y x hf y y +=+ Step 4 ),(1n n n n y x hf y y +=+ Euler 法解常微分方程算程序： function euler2(fun,y0,A,h) %fun--y' %y0---初值 %A----x 取值范围 %a----x 左区间端点值 %b----x 右区间端点值 %h----给定步长 x=min(A); b=max(A); y=y0; while x 指导教师： 年 月 日 改进Euelr 法解常微分方程 改进Euler 法解常微分方程算法： Step 1 分别取积分上限、积分下限、步长 Step 2 取一个以h 为步长，a ，b 分别为左右端点的矩阵 Step 3 （1）做显性Euler 预测),(1n n i i y x hf y y +=+ （2）将1+i y 带入)],(),([2 h 111+++++=i i i i i i y x f y x f y y Step 4计算h n n +=判断b n ≤是否成立，成立返回Step 3，否则继续进行Step 5 Step 5 )],(),([2 h 111+++++=i i i i i i y x f y x f y y 改进Euler 法解常微分方程算程序： function gaijineuler2(fun,y0,A,h) %fun--y' %y0---初值 %A----x 取值范围 %a----x 左区间端点值 %b----x 右区间端点值 %h----给定步长 a=min(A); b=max(A); x=a:h:b; y(1)=y0; for i=1:length(x)-1 w1=feval(fun,x(i),y(i)); y(i+1)=y(i)+h*w1; w2=feval(fun,x(i+1),y(i+1)); y(i+1)=y(i)+h*(w1+w2)/2; end x=x' 微分方程求解是系统仿真、数学模型实现以及很多工程问题求解的核心部分，应用MATLAB可以方便地对一阶常微分方程组进行求解，这里将对其基本方法进行介绍。值得注意的是，高阶微分方程组可以通过引进参变量化为一阶常微分方程组，也可以同样方便解决。 若有一个微分方程（组）的参变量为列向量,即,且它参变量随时间变化的微分方程可以有以下方程描述： 这里的f函数是一个列向量，即, i=1,2,3…,n，它可以是任意非线性函数。 则一般微分方程可以如此求解： [t,x]=ode45(f,timespan,x0) 对于刚性方程，即一些解变化缓慢，一些解变化剧烈，且两者相差较为悬殊的这种方程，通常调用ode15s而非o de45进行求解。 例1： 解：编写function或者用匿名函数 表达f=y-2*x/y即可； function dy=f(t,y) dy=y-2*t/y; end 命令： t=[0,1];%y0=1; [x,y]=ode45('f',t,1);%注意 这里的x相当于自变量t plot(x,y,x,sqrt(1+2*x)),legend('数值解','解析解'); 可见求解效果不错。 例2、 解：编写function function dx=f(t,x)%返回值是列向量 dx=[-x(2)-x(3); x(1)+0.2*x(2); 0.2+(x(1)-5.7)*x(3)]; end 命令： t=[0,100]; y0=[0 0 0]';%注意是列向量 [x,y]=ode45('f',t,y0); plot(x,y)； 例3、 这是一个二阶微分方程组，可以引进变量,由此ODE可以化成如下形式 可以采用和例2相同的方法求解： function dx=f(t,x) dx=[x(2); -(x(1)^2-1)*x(2)-x(1)]; End 第五讲 Matlab求解微分方程 教学目的：学会用MATLAB求简单微分方程的解析解、数值解，加深对微分方程概念和应用的理解；针对一些具体的问题，如追击问题，掌握利用软件求解微分方程的过程；了解微分方程模型解决问题思维方法及技巧；体会微分方程建摸的艺术性． 教学重点：利用机理分析建模微分方程模型，掌握追击问题的建模方法，掌握利用MATLAB求解数值解． 教学难点：利用机理分析建模微分方程模型，通过举例，结合图形以及与恰当的假设突破教学难点． 1微分方程相关函数（命令）及简介 因为没有一种算法可以有效地解决所有的 ODE 问题，为此，Matlab 提供了多种求解器 Solver，对于不同的ODE 问题，采用不同的Solver． 要特别的是：ode23、ode45 是极其常用的用来求解非刚性的标准形式的一阶常微分方程(组)的初值问题的解的 Matlab 的常用程序，其中： ode23 采用龙格-库塔2 阶算法，用3 阶公式作误差估计来调节步长，具有低等的精度． ode45 则采用龙格-库塔4 阶算法，用5 阶公式作误差估计来调节步长，具有中等的精度． 2 求解微分方程的一些例子 几个可以直接用 Matlab 求微分方程精确解的例子： 例1：求解微分方程 22x xe xy dx dy -=+，并加以验证． 求解本问题的Matlab 程序为： syms x y %line1 y=dsolve('Dy+2*x*y=x*exp(-x^2)','x') %line2 diff(y,x)+2 *x*y-x*exp(-x^2) %line3 simplify(diff(y,x)+2*x*y-x*exp(-x^2)) %line4 说明：(1) 行line1是用命令定义x,y 为符号变量．这里可以不写，但为确保正确性，建议写上； (2) 行line2是用命令求出的微分方程的解： 1/2*exp(-x^2)*x^2+exp(-x^2)*C1 (3) 行line3使用所求得的解．这里是将解代入原微分方程，结果应该为0，但这里给出： -x^3*exp(-x^2)-2*x*exp(-x^2)*C1+2*x*(1/2*exp(-x^2)*x^2+exp(-x^ Matlab 求解常微分方程边值问题的方法：bvp4c 函数 常微分方程的边值问题，即boundary value problems ，简称BVP 问题，是指表达形式为 (,)((),())0'=??=?y f x y g y a y b 或(,,)((),(),)0'=??=? y f x y p g y a y b p 的方程组（p 是未知参数），在MA TLAB 中使用积分器bvp4c 来求解。 [命令函数] bvp4c [调用格式] sol=bvp4c(odefun,bcfun,solinit,options,p1,p2,…) sol 为一结构体，sol.x 、sol.y 、sol.yp 分别是所选择的网格点及其对应的y(x)与y'(x)数值; bvp4c 为带边值条件常微分方程积分器的函数命令；odefun 为描述微分方程组的函数文件；bcfun 为计算边界条件g(f(a),f(b),p)=0的函数文件；solinit 为一结构体，solinit.x 与solinit.y 分别是初始网格的有序节点与初始估计值，边界值条件分别对应a=solinit.x(l)和b=solinit.x(end)； options 为bvpset 命令设定的可选函数，可采用系统默认值；p1, p2…为未知参数。 例 求常微分方程0''+=y y 在(0)2=y 与(4)2=-y 时的数值解。 [解题过程] 仍使用常用方法改变方程的形式： 令1=y y ，21'=y y ，则原方程等价于标准形式的方程组1221 ?'=??'=-??y y y y ； 将其写为函数文件twoode.m ； 同时写出边界条件函数对应文件twobc.m ； 分别使用结构solinit 和命令bvp4c 确定y-x 的关系； 作出y-x 的关系曲线图。 [算例代码] solinit =bvpinit(linspace(0,4,5),[1 0]); % linspace(0,4,5)为初始网格，[1,0]为初始估计值 sol=bvp4c(@twoode,@twobc,solinit); % twoode 与twobc 分别为微分方程与边界条件的函数，solinit 为结构 x=linspace(0,4); %确定x 范围 y=deval(sol,x); %确定y 范围 plot(x,y(1,:)); %画出y-x 的图形 %定义twoode 函数（下述代码另存为工作目录下的twoode.m 文件） function dydx= twoode(x,y) %微分方程函数的定义 dydx =[y(2) -abs(y(1))]; %定义twobc 函数（下述代码另存为工作目录下的twobc.m 文件） function res= twobc(ya,yb); %边界条件函数的定义 res=[ya(1);yb(1)+2];Matlab求解微分方程(组)及偏微分方程(组)
Matlab解微分方程(ODE+PDE)
matlabeuler法解常微分方程
常微分方程组的MATLAB求解范例
第五讲matlab求解微分方程
Matlab求解常微分方程边值问题的方法
相关主题
文本预览