实验二微分方程与差分方程模型Matlab求解
matlab差分法求解微分方程
一、概述微分方程是自然科学和工程技术中常见的数学模型，它描述了连续系统的变化规律。
在实际应用中，求解微分方程是一项重要且复杂的工作。
而matlab是一种常用的科学计算软件，它提供了丰富的数学函数和工具，能够辅助工程师和科学家在求解微分方程方面取得良好的效果。
二、matlab差分法求解微分方程的基本原理差分法是一种常见的数值求解微分方程的方法。
它基于微分的定义，将微分方程中的微分运算用差分逼近来进行计算。
在matlab中，可以利用内置的数学函数和工具，通过差分法求解微分方程，得到数值解或者近似解。
三、matlab中使用差分法求解常微分方程的步骤1. 确定微分方程的类型和边界条件需要明确所要求解的微分方程是什么类型的，以及其所对应的边界条件是什么。
这对于后续的数值求解过程非常重要。
在matlab中，可以利用符号变量和函数来表示微分方程和边界条件。
2. 将微分方程离散化接下来，需要将微分方程进行离散化处理，将微分方程中的微分运算用差分逼近来进行计算。
这一步需要根据微分方程的具体形式和求解精度选择合适的差分方法，常见的有前向差分、后向差分和中心差分等方法。
3. 构建代数方程组将离散化后的微分方程转化为代数方程组。
这一步需要根据微分方程的离散化表达式和边界条件，利用matlab的矩阵和向量运算功能，构建代数方程组。
4. 求解代数方程组利用matlab的求解函数，求解构建得到的代数方程组，得到微分方程的数值解或者近似解。
在求解过程中，需要注意数值稳定性和收敛性，以及选择合适的数值积分方法和迭代算法。
四、实例：使用matlab差分法求解一阶常微分方程为了更好地理解matlab中使用差分法求解微分方程的过程，以下将通过一个具体的实例来演示。
假设要求解如下的一阶常微分方程：dy/dx = -2x + 1, y(0) = 11. 确定微分方程的类型和边界条件根据给定的方程，可以确定它是一阶常微分方程，且给定了初始条件y(0) = 1。
实验二微分方程与差分方程模型Matlab求解
实验二： 微分方程与差分方程模型Matlab 求解一、实验目的[1] 掌握解析、数值解法，并学会用图形观察解的形态和进展解的定性分析；[2] 熟悉MATLAB 软件关于微分方程求解的各种命令； [3] 通过范例学习建立微分方程方面的数学模型以及求解全过程；[4] 熟悉离散 Logistic 模型的求解与混沌的产生过程。
二、实验原理1. 微分方程模型与MATLAB 求解 解析解用MATLAB 命令dsolve(‘eqn1’,’eqn2’, ...) 求常微分方程〔组〕的解析解。
其中‘eqni'表示第i 个微分方程，Dny 表示y 的n 阶导数,默认的自变量为t 。
〔1〕 微分方程例1 求解一阶微分方程 21y dxdy+= (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()04(/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 。
MATLAB差异方程与微分方程求解技巧
MATLAB差异方程与微分方程求解技巧差异方程和微分方程是数学中重要的概念和工具，它们在各个领域的建模和分析中发挥着重要作用。
而MATLAB作为一款强大的数学软件，提供了丰富的工具和函数来求解差异方程和微分方程。
本文将介绍MATLAB中差异方程和微分方程的求解技巧，并提供一些实际案例来加深理解。
一、差异方程的求解技巧差异方程是描述离散域系统的数学模型，通常用递归关系来表达。
MATLAB 提供了多种方法来求解差异方程，其中最常用的是通过递推关系进行迭代。
1. 递推法递推法是通过迭代计算差异方程中的每一项来求解整个方程。
首先，需要定义差异方程的初始条件和递推关系。
然后，可以使用循环结构来进行迭代计算，直到达到所需精度或迭代次数。
假设我们要求解以下差异方程：y[n] = a * y[n-1] + b * y[n-2]其中，a和b为常数，y[n]为求解的项，y[n-1]和y[n-2]为已知的前两项。
在MATLAB中，可以使用for循环或while循环来实现递推法求解差异方程。
以下是使用for循环的实例代码：``` MATLABn = 1:10; % 定义计算的范围y = zeros(size(n)); % 初始化y的空间y(1) = y0; % 设定初始条件y(2) = y1; % 设定初始条件for i = 3:length(n)y(i) = a * y(i-1) + b * y(i-2); % 递推计算end```2. 齐次差异方程和非齐次差异方程的求解在求解差异方程时，需要区分齐次差异方程和非齐次差异方程。
对于齐次差异方程，它的非零解为零解；对于非齐次差异方程，它的非零解可以通过叠加齐次解和特解来得到。
MATLAB中，可以使用dsolve函数来求解差异方程。
以下是求解一阶齐次差异方程的实例代码：``` MATLABsyms y(t); % 定义符号变量eqn = diff(y, t) == a * y; % 定义差异方程cond = y(0) == y0; % 定义初始条件ySol(t) = dsolve(eqn, cond); % 求解差异方程```二、微分方程的求解技巧微分方程是描述连续域系统的数学模型，通常用导数关系来表达。
实验二_基于Matlab的微分方程数值解法
实验二微分方程数值解法一．实验原理及实验内容：对微分方程描述的控制系统，利用欧拉法、二阶龙格-库塔法、四阶龙格-库塔法分别编写M文件，进行数值计算和作图。
1．分别用欧拉法、二阶龙格-库塔法、四阶龙格-库塔法求下面系统的输出响应y(t)在0≤t≤1上，h=0.1时的数值解。
'2,(0)1=-=y y y要求保留4位小数，并将三种方法的结果与真解2=进行比较。
()ty t e-2.若为如2y y y==何编程计算？',(0)1二．实验仪器：计算机Matlab软件三．实验数据记录：程序一：disp('欧拉算法');y=1;h=0.1;for i=0:0.1:1disp(y);y=y+h*(-2*y);enddisp('欧拉算法');ydisp('精确解');yy=exp(-2*t)h=0.1;disp('函数的2阶数值解为');disp('y=');y=1;for t=0:h:1;disp(y);k1=-2*y;k2=-2*(y+k1*h);y(i+1)=y(i)+(k1+k2)*h*1/2;endh=0.1;disp('函数的4阶数值解为');disp('y=');y=1;for t=0:h:1;disp(y);k1=-2*y;k2=-2*(y+k1*h*1/2);k3=-2*(y+k2*h*1/2);k4=-2*(y+k3*h);y=y+h*1/6*(k1+2*k2+2*k3+k4); end>>程序2：t=0:0.1:1;n=length(t);y(1)=1;h=0.1;for i=1:n-1y(i+1)=y(i)+h*(y(i)*y(i)); enddisp('欧拉算法');ydisp('精确解');yy=exp(-2*t)h=0.1;disp('函数的2阶数值解为');disp('y=');y=1;for t=0:h:1;disp(y);k1=y*y;k2=(y+k1*h)^2;y=y+(k1+k2)*h*1/2;endh=0.1;disp('函数的4阶数值解为');disp('y=');y=1;disp(y);k1=y*y;k2=(y+k1*h*1/2)^2;k3=(y+k2*h*1/2)^2;k4=(y+k3*h)^2;y=y+h*1/6*(k1+2*k2+2*k3+k4); end。
数学建模实验二：微分方程与差分方程模型Matlab求解
实验二：微分方程与差分方程模型Matlab 求解专业年级： 2014级信息与计算科学1班 姓名： 黄志锐 学号：201430120110一、实验目的[1] 掌握解析、数值解法，并学会用图形观察解的形态和进行解的定性分析； [2] 熟悉MATLAB 软件关于微分方程求解的各种命令；[3] 通过范例学习建立微分方程方面的数学模型以及求解全过程； [4] 熟悉离散 Logistic 模型的求解与混沌的产生过程。
二、实验内容1．求微分方程的解析解, 并画出图形，解：使用MATLAB 编程计算得出上述微分方程的解析解为：y =3e x −2x −2其解析解图示如图1所示：图1 解析解图示=+2,(0)1,01y y x y x '=<<MATLAB代码运行结果截图如下所示：2．求微分方程的数值解, 并画出图形，解：使用MATLAB编程计算上述微分方程的数值解，并作出其数值解的图示如图2所示：图2 数值解图示cos0,(0)1,(0)0y y x y y'''+===3．两种相似的群体之间为了争夺有限的同一种食物来源和生活空间而进行生存竞争时，往往是竞争力较弱的种群灭亡，而竞争力较强的种群达到环境容许的最大数量。
假设有甲、乙两个生物种群，当它们各自生存于一个自然环境中，均服从 Logistic 规律。
（1）是两个种群的数量； （2）是它们的固有增长率； （3）是它们的最大容量；（4）为种群乙(甲)占据甲(乙)的位置的数量，并且建立模型如下：)(),(21t x t x 21,r r 21,n n )(12m m .1122;x m x m βα==计算, 画出图形及相轨线图。
解释其解变化过程。
2）121,r r ==设12100,n n ==102010x x ==,=1.5，=0.7，计算, 画出图形及相轨线图。
解释其解变化过程。
解：（1）使用MATLAB 编程实现上述模型，并输入相关参数后，计算得出, 并画出图形及相轨线图如下所示：图3 数值解图示)(),(21t x t x αβ)(),(21t x t x )(),(21t x tx图4 相轨线图详细MATLAB代码如下：由图3、图4可以看出，甲、乙两个生物种群以几乎一致的趋势不断增长，直到达到一个相对稳定的数量。
用MATLAB求解微分方程及微分方程组
例 3 求微分方程组的通解. dx dt 2 x 3 y 3 z dy 4 x 5 y 3z dt dz 4 x 4 y 2 z dt
任取k1、k2的一组初始值：k0=[2,1];
输入命令： k=lsqcurvefit('curvefun1',k0,t,c) 运行结果为： k =[ 1.3240 作图表示求解结果： t1=0:0.1:18; f=curvefun1(k,t1); plot(t,c,'ko',t1,f,'r-')
90 80 70 60 50 40 30 20 10 0
0.2573]
0
2
46Leabharlann 81012
14
16
18
模型2：慢速饮酒时，体液中酒精含量的变化率
dx k2 x a dt x(0) 0
其中
M a T
M为饮酒的总量，T为饮酒的时间
则有；
a x (1 e k 2 t ) k2
5 5 ) 处时被导弹击中. 当 x 1时 y ，即当乙舰航行到点 (1, 24 24 y 5 被击中时间为: t . 若 v0=1, 则在 t=0.21 处被击中. v0 24v0
轨迹图如下
例： 饮酒模型
模型1：快速饮酒后，胃中酒精含量的变化率
dy k1 y dt y (0) M
5 5 ) 处时被导弹击中. 当 x 1时 y ，即当乙舰航行到点 (1, 24 24 y 5 被击中时间为 : t . 若 v0=1, 则在 t=0.21 处被击中. v0 24v0
实验二MATLAB 求微分方程的解
实验二 微分方程求解一、问题背景与实验目的实际应用问题通过数学建模所归纳而得到的方程，绝大多数都是微分方程，真正能得到代数方程的机会很少．另一方面，能够求解的微分方程也是十分有限的，特别是高阶方程和偏微分方程（组）．这就要求我们必须研究微分方程（组）的解法，既要研究微分方程（组）的解析解法（精确解），更要研究微分方程（组）的数值解法（近似解）．对微分方程（组）的解析解法(精确解)，Matlab 有专门的函数可以用，本实验将作一定的介绍．本实验将主要研究微分方程(组)的数值解法（近似解），重点介绍 Euler 折线法．二、相关函数（命令）及简介1．dsolve('equ1','equ2',…)：Matlab 求微分方程的解析解．equ1、equ2、…为方程（或条件）．写方程（或条件）时用 Dy 表示y 关于自变量的一阶导数，用用 D2y 表示 y 关于自变量的二阶导数，依此类推．2．simplify(s)：对表达式 s 使用 maple 的化简规则进行化简． 例如： syms xsimplify(sin(x)^2 + cos(x)^2) ans=13．[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 = combine4．[T,Y] = solver(odefun,tspan,y 0) 求微分方程的数值解． 说明：(1) 其中的 solver 为命令 ode45、ode23、ode113、ode15s 、ode23s 、ode23t 、ode23tb 之一．(2) odefun是显式常微分方程：⎪⎩⎪⎨⎧==00)(),(yt y y t f dt dy(3) 在积分区间 tspan =],[0f t t 上，从0t 到f t ，用初始条件0y 求解．(4) 要获得问题在其他指定时间点 ,210,,t t t 上的解，则令 tspan =],,,[,210f t t t t （要求是单调的）．(5) 因为没有一种算法可以有效地解决所有的 ODE 问题，为此，Matlab 提供了多种求解器 Solver ，对于不同的ODE 问题，采用不同的Solver ．(6) 要特别的是：ode23、ode45 是极其常用的用来求解非刚性的标准形式的一阶常微分方程(组)的初值问题的解的 Matlab 的常用程序，其中：ode23 采用龙格-库塔2 阶算法，用3 阶公式作误差估计来调节步长，具有低等的精度．ode45 则采用龙格-库塔4 阶算法，用5 阶公式作误差估计来调节步长，具有中等的精度．5．ezplot(x,y ,[tmin,tmax])：符号函数的作图命令．x,y 为关于参数t 的符号函数，[tmin,tmax] 为 t 的取值范围．6．inline()：建立一个内联函数．格式：inline('expr', 'var1', 'var2',…) ，注意括号里的表达式要加引号．例：Q = dblquad(inline('y*sin(x)'), pi, 2*pi, 0, pi)三、实验内容1. 几个可以直接用 Matlab 求微分方程精确解的例子： 例1：求解微分方程22xxexy dxdy -=+，并加以验证．求解本问题的Matlab 程序为：syms x y %line1 y=dsolve('Dy+2*x*y=x*exp(-x^2)','x') %line2diff(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^2)*C1)(4) 行line4 用 simplify() 函数对上式进行化简，结果为 0， 表明)(x y y =的确是微分方程的解．例2：求微分方程0'=-+x e y xy 在初始条件e y 2)1(=下的特解，并画出解函数的图形．求解本问题的 Matlab 程序为： syms x yy=dsolve('x*Dy+y-exp(x)=0','y(1)=2*exp(1)','x') ezplot(y)微分方程的特解为：y=1/x*exp(x)+1/x* exp (1) (Matlab 格式)，即xe e y x+=，解函数的图形如图 1：图1例3：求微分方程组⎪⎪⎩⎪⎪⎨⎧=--=++035y x dt dy e y x dtdx t在初始条件0|,1|00====t t y x 下的特解，并画出解函数的图形．求解本问题的 Matlab 程序为： syms x y t[x,y]=dsolve('Dx+5*x+y=exp(t)','Dy-x-3*y=0','x(0)=1','y(0)=0','t') simple(x); simple(y);ezplot(x,y ,[0,1.3]);axis auto微分方程的特解(式子特别长)以及解函数的图形均略． 2. 用ode23、ode45等求解非刚性的标准形式的一阶常微分方程(组)的初值问题的数值解(近似解)．例4：求解微分方程初值问题⎪⎩⎪⎨⎧=++-=1)0(2222y x x y dxdy 的数值解，求解范围为区间[0, 0.5]．fun=inline('-2*y+2*x^2+2*x','x','y'); [x,y]=ode23(fun,[0,0.5],1); x'; y';plot(x,y ,'o-') >> x' ans =0.0000 0.0400 0.0900 0.1400 0.1900 0.2400 0.2900 0.3400 0.3900 0.4400 0.4900 0.5000 >> y' ans =1.0000 0.9247 0.8434 0.7754 0.7199 0.6764 0.6440 0.6222 0.6105 0.6084 0.6154 0.6179 图形结果为图 2．图2例 5：求解描述振荡器的经典的 V er der Pol 微分方程.7,0)0(',1)0(,0)1(222====+--μμy y y dtdy y dty d分析：令,,121dtdx x y x ==则.)1(,1221221x x x dtdx x dtdx --==μ先编写函数文件verderpol.m ： function xprime = verderpol(t,x) global mu;xprime = [x(2);mu*(1-x(1)^2)*x(2)-x(1)]; 再编写命令文件vdp1.m ： global mu; mu = 7; y0=[1;0][t,x] = ode45('verderpol',[0,40],y0); x1=x(:,1);x2=x(:,2); plot(t,x1)图形结果为图3．图33. 用 Euler 折线法求解前面讲到过，能够求解的微分方程也是十分有限的．下面介绍用 Euler 折线法求微分方程的数值解(近似解)的方法．Euler 折线法求解的基本思想是将微分方程初值问题⎪⎩⎪⎨⎧==00)(),,(yx y y x f dx dy化成一个代数方程，即差分方程，主要步骤是用差商hx y h x y )()(-+替代微商dxdy ,于是：⎪⎩⎪⎨⎧==-+)()),(,()()(00x y y x y x f h x y h x y k k k k 记)(,1k k k k x y y h x x =+=+，从而)(1h x y y k k +=+，则有1,,2,1,0).,(,),(1100-=⎪⎩⎪⎨⎧+=+==++n k y x hf y yh x x x y y k k k k k k 例 6：用 Euler 折线法求解微分方程初值问题⎪⎩⎪⎨⎧=+=1)0(,22y y x y dxdy 的数值解(步长h 取0.4)，求解范围为区间[0,2]．解：本问题的差分方程为1,,2,1,0).2),( ),(,,4.0,1,021100-=⎪⎪⎪⎩⎪⎪⎪⎨⎧+=+=+====++n k y x y y x f y x hf y y h x x h y x k k k k k k （其中： 相应的Matlab 程序见附录 1． 数据结果为：0 1.0000 0.4000 1.4000 0.8000 2.1233 1.2000 3.1145 1.6000 4.4593 2.0000 6.3074图形结果见图4：图4特别说明：本问题可进一步利用四阶 Runge-Kutta 法求解，读者可将两个结果在一个图中显示，并和精确值比较，看看哪个更“精确”？（相应的 Matlab 程序参见附录 2）．四、自己动手1. 求微分方程0sin 2')1(2=-+-x xy y x 的通解．2. 求微分方程x e y y y x sin 5'2''=+-的通解．3. 求微分方程组⎪⎪⎩⎪⎪⎨⎧=-+=++00y x dtdy y x dtdx在初始条件0|,1|00====t t y x 下的特解，并画出解函数()y f x =的图形． 4. 分别用 ode23、ode45 求上述第 3 题中的微分方程初值问题的数值解(近似解)，求解区间为[0,2]t ∈．利用画图来比较两种求解器之间的差异．5. 用 Euler 折线法求解微分方程初值问题⎪⎩⎪⎨⎧=-=1)0(,12'32y y xy y 的数值解(步长h 取0.1)，求解范围为区间[0,2]．6. 用四阶 Runge-Kutta 法求解微分方程初值问题⎩⎨⎧=-=1)0(,cos 'y x e y y x 的数值解(步长h 取0.1)，求解范围为区间[0,3]．四阶 Runge-Kutta 法的迭代公式为(Euler 折线法实为一阶 Runge-Kutta 法)：1,,2,1,0),()2,2()2,2(),()22(6,),(342312143211100-=⎪⎪⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎪⎪⎨⎧++=++=++==++++=+==++n k hL y h x f L L h y h x f L L h y h x f L y x f L L L L L hy y h x x x y y k k k k k k k k k k k k 相应的 Matlab 程序参见附录 2．试用该方法求解第5题中的初值问题． 7. 用 ode45 方法求上述第 6 题的常微分方程初值问题的数值解(近似解)，从而利用画图来比较两者间的差异．五、附录附录 1：(fulu1.m)clearf=sym('y+2*x/y^2'); a=0; b=2; h=0.4;n=(b-a)/h+1; x=0; y=1;szj=[x,y]; for i=1:n-1y=y+h*subs(f,{'x','y'},{x,y}); x=x+h;szj=[szj;x,y]; end szjplot(szj(:,1),szj(:,2))附录 2：(fulu2.m)clearf=sym('y-exp(x)*cos(x)'); a=0; b=3; h=0.1;n=(b-a)/h+1; x=0; y=1;szj=[x,y];for i=1:n-1l1=subs(f,{'x','y'},{x,y});l2=subs(f,{'x','y'},{x+h/2,y+l1*h/2});l3=subs(f,{'x','y'},{x+h/2,y+l2*h/2});l4=subs(f,{'x','y'},{x+h,y+l3*h});y=y+h*(l1+2*l2+2*l3+l4)/6;x=x+h;szj=[szj;x,y];endszjplot(szj(:,1),szj(:,2))。
matlab差分法解微分方程
matlab差分法解微分方程在MATLAB中，差分法是一种常用的数值方法，用于解决微分方程。
差分法的基本思想是将微分方程中的导数用离散的差分近似表示，然后通过迭代计算得到方程的数值解。
下面我将从多个角度来解释如何使用差分法在MATLAB中解微分方程。
1. 离散化，首先，我们需要将微分方程离散化，将自变量和因变量分成若干个离散的点。
例如，可以选择一个均匀的网格，将自变量的取值离散化为一系列的点。
这样，微分方程中的导数可以用差分近似来表示。
2. 差分近似，使用差分近似来代替微分方程中的导数。
最常见的差分近似方法是中心差分法。
对于一阶导数，可以使用中心差分公式，f'(x) ≈ (f(x+h) f(x-h)) / (2h)，其中h是离散化步长。
对于二阶导数，可以使用中心差分公式，f''(x) ≈ (f(x+h) 2f(x) + f(x-h)) / (h^2)。
根据微分方程的类型和边界条件，选择适当的差分近似方法。
3. 矩阵表示，将差分近似后的微分方程转化为矩阵形式。
通过将微分方程中的各项离散化，可以得到一个线性方程组。
这个方程组可以用矩阵表示，其中未知量是离散化后的因变量。
4. 数值求解，使用MATLAB中的线性代数求解函数，例如backslash运算符（\）或者LU分解等，求解得到线性方程组的数值解。
这个数值解就是微分方程的近似解。
需要注意的是，差分法是一种数值方法，所得到的解是近似解，精确度受离散化步长的影响。
通常情况下，可以通过减小离散化步长来提高数值解的精确度。
此外，对于某些特殊类型的微分方程，可能需要采用更高级的差分方法，如龙格-库塔法（Runge-Kutta method）或有限元方法（Finite Element Method）等。
综上所述，差分法是一种常用的数值方法，可以在MATLAB中用于解决微分方程。
通过离散化、差分近似、矩阵表示和数值求解等步骤，可以得到微分方程的数值解。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
实验二: 微分方程与差分方程模型Matlab 求解一、实验目的[1] 掌握解析、数值解法,并学会用图形观察解的形态与进行解的定性分析; [2] 熟悉MATLAB 软件关于微分方程求解的各种命令;[3] 通过范例学习建立微分方程方面的数学模型以及求解全过程; [4] 熟悉离散 Logistic 模型的求解与混沌的产生过程。
二、实验原理1、 微分方程模型与MATLAB 求解解析解用MATLAB 命令dsolve(‘eqn1’,’eqn2’, 、、、) 求常微分方程(组)的解析解。
其中‘eqni'表示第i 个微分方程,Dny 表示y 的n 阶导数,默认的自变量为t 。
(1) 微分方程例1 求解一阶微分方程 21y dxdy+= (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()04(/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 子目录,也可以保存在自定义文件夹,这时注意要增加搜索路径( 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。
+uu''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、符号说明:①x1(t), x2(t)分别就是食饵、捕食者(鲨鱼)在t时刻的数量;②r1, r2就是食饵、捕食者的固有增长率;③λ1就是捕食者掠取食饵的能力,λ2就是食饵对捕食者的供养能力;2、基本假设:① 捕食者的存在使食饵的增长率降低,假设降低的程度与捕食者数量成正比,即)(21111x r x dtdx λ-= ② 食饵对捕食者的数量x 2起到增长的作用, 其程度与食饵数量x 1成正比,即)(12222x r x dtdx λ+-= 综合以上①与②,得到如下模型:模型一:不考虑人工捕获的情况该模型反映了在没有人工捕获的自然环境中食饵与捕食者之间的制约关系,没有考虑食饵与捕食者自身的阻滞作用,就是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⎪⎪⎩⎪⎪⎨⎧+-=-=)()(1222221111x r x dtdx x r x dt dx λλ2)0(,25)0()02.05.0()1.01(21122211==⎪⎪⎩⎪⎪⎨⎧+-=-=x x x x dtdx x x dt dx从图中可以瞧出它们的数量都呈现出周期性,而且鲨鱼数量的高峰期稍滞后于食饵数量的高峰期。
画出相轨迹图:假设人工捕获能力系数为e,相当于食饵的自然增长率由r 1 降为r 1-e,捕食者的死亡率由r 2 增为 r 2+e,因此模型一修改为:设战前捕获能力系数e=0、3, 战争中降为e=0、1, 其它参数与模型(一)的参数相同。
观察结果会如何变化？⎪⎩⎪⎨⎧++-=--=])([])[(1222221111x e r x dtdxx 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);plot(t1,p1,'-',t2,p2,'*')11222112(0.70.1)(0.80.02)(0)25,(0)2dx x x dt dx x x dtx x ⎧=-⎪⎪⎪=-+⎨⎪==⎪⎪⎩⎪⎪⎪⎩⎪⎪⎪⎨⎧==+-=-=2)0(,25)0()02.06.0()1.09.0(21122211y y y y dtdy y y dt dy )()()()(;)()()()(21222121t y t y t y t p t x t x t x t p +=+=0510150.10.20.30.40.50.60.70.8图中‘*’曲线为战争中鲨鱼所占比例。