数学建模实验二:微分方程模型Matlab求解与分析
用matlab对微分方程求解实验报告.
o 《高等数学》上机作业(三一、上机目的1、学会用 M a t l a b 求简单微分方程的解析解。
2、学会用 M a t l a b 求微分方程的数值解。
二、上机内容1、求简单微分方程的解析解.2、求微分方程的数值解.3、数学建模实例.4、上机作业. 三、上机作业1. 求微分方程:在初值条件下的特解,并画出解函数的图形. 命令>> y =d s o l v e ('x *D y +y -e x p (x =0','y (1=2*e x p (1','x ' 运行结果:y = 1/x *e x p (x +1/x *e x p (1'xxy y e +-=12(y e =函数图象:2. 求微分方程的特解.22450(00,'(110d y dyy dx dx y y ?+-=???==?命令>> y=dsolve('D2y+4*Dy-5*y=0','y(0=0,Dy(1=10','x' 运行结果:y=10/(exp(1+5*exp(-5*exp(x-10/(exp(1+5*exp(-5*exp(-5*x3. 鱼雷追击问题一敌舰在某海域内沿着正北方向航行时,我方战舰恰好位于敌舰的正西方向 1 公里处.我舰向敌舰发射制导鱼雷,敌舰速度为0.42 公里/分,鱼雷速度为敌舰速度的2倍。
试问敌舰航行多远时将被击中?M文件x0=0; xf=0.9999999999999; [x,y]=ode15s('eq1',[x0 xf],[0 0]; plot(x,y(:,1,'b.'hold on;y=0:0.1:1;plot(1,y, '*'运行结果图像:结论:大概在y=0.67处击中敌方舰艇!(选做一个慢跑者在平面上沿椭圆以恒定的速率v=1跑步,设椭圆方程为:x=10+20cost, y=20+5sint. 突然有一只狗攻击他. 这只狗从原点出发,以恒定速率w跑向慢跑者,狗的运动方向始终指向慢跑者.分别求出w=20,w=5时狗的运动轨迹.W=20M文件代码function dy=eq3(t,ydy=zeros(2,1;dy(1=20*(10+20*cos(t-y(1/sqrt((10+20*cos(t-y(1^2+(2 0+15*sin(t-y(2^2;dy(2=20*(20+15*sin(t-y(2/sqrt((10+20*cos(t-y(1^2+(2 0+15*sin(t-y(2^2;运行命令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 onplot(y(:,1,y(:,2,'r*'运行结果:利用二分法更改tf tf=5时tf=2.5时tf=3.15时:所以在t=3.15时刻恰好追上!W=5M文件代码function dy=eq4(t,ydy=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;[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 onplot(y(:,1,y(:,2,'*' 运行结果更改tf=20运行结果Tf=40 11所以永远追不上!四、上机心得体会高等数学是工科学生的主干科目,它应用于生产生活的方方面面,通过建模,计算可以求出实际问题的最优化问题!因此我们需要掌握建模和利用专业软件处理实际问题的能力! 12。
微分方程模型matlab
微分方程模型matlab微分方程是自然科学及工程领域中常用的数学工具,它描述了系统中的变化率与其当前状态之间的关系。
而Matlab是一个广泛使用的科学计算软件,它具有强大的数值计算能力和可视化功能。
在本文中,我们将讨论如何使用Matlab来求解微分方程模型,并给出一个实例。
一、Matlab中的微分方程求解函数Matlab中有几个可用于求解微分方程的函数,包括ode45、ode23、ode113、ode15s等。
其中前两者是比较常用的。
这些函数都采用了一些数值方法来近似求解微分方程,得到数值解。
这些求解函数都采用输入函数来表示微分方程。
二、小例子:一个单自由度振动模型作为一个简单的微分方程模型,我们来考虑一个单自由度振动模型。
它可以用下面的微分方程来表示:$$m\ddot{x}+c\dot{x}+kx=F(t)$$其中,$m$是质量,$c$是阻尼系数,$k$是弹簧刚度,$F(t)$是作用在系统上的外力,$x(t)$是位移,$\dot{x}(t)$是速度,$\ddot{x}(t)$是加速度。
我们假设质量$m$为1kg,阻尼系数$c$为0.2N-s/m,弹簧刚度$k$为100N/m,外力$F(t)$为3N。
首先,我们将这些参数定义为变量。
```matlabm = 1;c = 0.2;k = 100;F = 3;```然后,我们将微分方程表示为Matlab可以识别的形式。
我们定义一个名为odefun的函数,在此函数中计算微分方程的右侧。
```matlabfunction dxdt = odefun(t, x)dxdt = [x(2); (F - c*x(2) - k*x(1))/m];end```接着,我们要设置初始状态。
我们初始化$x_0$为0,$\dot{x}_0$为0.5m/s。
```matlabx0 = [0; 0.5];```我们使用ode45函数来解微分方程。
ode45函数需要我们提供微分方程函数和初始条件,然后它会自动计算微分方程的解。
使用Matlab进行微分方程求解的方法
使用Matlab进行微分方程求解的方法引言微分方程是数学中非常重要的一部分,广泛应用于物理、经济、工程等领域。
对于大部分微分方程的解析解往往难以求得,而数值解法则成为了一种常用的解决手段。
Matlab作为一种强大的科学计算软件,也提供了丰富的工具和函数用于求解微分方程,本文将介绍一些常见的使用Matlab进行微分方程求解的方法。
一、数值求解方法1. 欧拉方法欧拉方法是最简单的一种数值求解微分方程的方法,它将微分方程的微分项用差分的方式进行近似。
具体的公式为:y(n+1) = y(n) + hf(x(n), y(n))其中,y(n)表示近似解在第n个点的值,h为步长,f(x, y)为微分方程的右端项。
在Matlab中使用欧拉方法进行求解可以使用ode113函数,通过设定不同的步长,可以得到不同精度的数值解。
2. 中点法中点法是较为精确的一种数值求解微分方程的方法,它的计算公式为:k1 = hf(x(n), y(n))k2 = hf(x(n) + h/2, y(n) + k1/2)y(n+1) = y(n) + k2中点法通过计算两个斜率的平均值来得到下一个点的值,相较于欧拉方法,中点法能提供更精确的数值解。
3. 4阶龙格库塔法龙格库塔法是一类高阶数值求解微分方程的方法,其中4阶龙格库塔法是最常用的一种。
它的计算公式为:k1 = hf(x(n), y(n))k2 = hf(x(n) + h/2, y(n) + k1/2)k3 = hf(x(n) + h/2, y(n) + k2/2)k4 = hf(x(n) + h, y(n) + k3)y(n+1) = y(n) + (k1 + 2k2 + 2k3 + k4)/64阶龙格库塔法通过计算多个斜率的加权平均值来得到下一个点的值,相较于欧拉方法和中点法,它的精度更高。
二、Matlab函数和工具除了可以使用以上的数值方法进行微分方程求解之外,Matlab还提供了一些相关的函数和工具,方便用户进行微分方程的建模和求解。
实验二_基于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微分方程模型
matlab微分方程模型Matlab微分方程模型是一种基于Matlab软件的数学建模方法,用于解决微分方程相关的问题。
微分方程是描述物理、工程和数学问题的重要工具,通过建立微分方程模型,可以对各种现象进行定量分析和预测。
在Matlab中,可以使用ode45函数求解常微分方程(ODE)或者ode15s函数求解刚性ODE。
这些函数可以通过数值方法近似求解微分方程的解析解,从而得到问题的数值解。
具体来说,可以通过在Matlab中定义微分方程的右侧函数,然后使用相应的ode函数进行求解。
例如,考虑一个简单的一阶线性微分方程模型:dy/dx = -ky,其中k为常数。
我们可以通过在Matlab中定义这个微分方程的右侧函数,并使用ode45函数求解。
具体步骤如下:1. 在Matlab中定义微分方程的右侧函数:function dydx = myODE(x,y)k = 0.1; % 设定常数k的值dydx = -k*y;end2. 使用ode45函数求解微分方程:xspan = [0 10]; % 设定求解区间y0 = 1; % 设定初始条件[x,y] = ode45(@myODE, xspan, y0);3. 绘制得到的数值解:plot(x,y);xlabel('x');ylabel('y');title('Solution of dy/dx = -ky');通过以上步骤,我们可以得到微分方程dy/dx = -ky的数值解,并绘制出解的图像。
这个简单的例子展示了如何使用Matlab微分方程模型求解微分方程。
除了一阶线性微分方程,Matlab微分方程模型还可以用于解决更复杂的微分方程问题,包括高阶线性微分方程、非线性微分方程、偏微分方程等。
通过定义相应的微分方程函数和合适的求解方法,可以在Matlab中进行数值求解。
此外,Matlab还提供了丰富的绘图和分析工具,可以对微分方程的解进行可视化和进一步分析。
实验二MATLAB数值计算常微分方程(组)的求解
实验⼆MATLAB数值计算常微分⽅程(组)的求解实验⼆ MATLAB 数值计算:常微分⽅程(组)的求解⼀、实验⽬的在物理学和⼯程技术上,很多问题都可以⽤⼀个或⼀组常微分⽅程来描述,因此要解决相应的实际问题往往需要⾸先求解对应的微分⽅程。
在⼤多数情况下这些微分⽅程通常是⾮线性的或者是超越⽅程(⽐如范德堡⽅程,波导本征值⽅程等),因此往往需要使⽤计算机数值求解。
MATLAB 作为⼀种强⼤的科学计算语⾔,其在数值计算和数据的可视化⽅⾯具有⽆以伦⽐的优势。
在解决常微分⽅程问题上,MATLAB 就提供了多种可适⽤于不同场合(如刚性和⾮刚性问题)下的求解器(Solver),例如ode45,ode15s ,ode23,ode23s 等等。
本次实验将以范德堡⽅程的计算和地球卫星的运⾏轨道的仿真为例,练习使⽤MATLAB 的常微分⽅程求解器,以期达到如下⼏个⽬的:1. 熟悉常微分⽅程的求解⽅法,了解状态⽅程的概念;2. 能熟练使⽤dsolve 函数解析求解常微分⽅程;3. 能熟练运⽤ode45、ode15s 求解器分别数值求解⾮刚性和刚性常微分⽅程;4. 学习⽤求解器来绘制相图的⽅法。
⼆、实验的预备知识1.微分⽅程的概念未知的函数以及它的某些阶的导数连同⾃变量都由⼀已知⽅程联系在⼀起的⽅程称为微分⽅程。
如果未知函数是⼀元函数,称为常微分⽅程(Ordinary differential equations ,简称odes )。
n 阶常微分⽅程的⼀般形式(隐式)为:0),,",',,()(=n y y y y t F (1)其中t 为⾃变量。
如果未知函数是多元函数,成为偏微分⽅程。
联系⼀些未知函数的⼀组微分⽅程组称为微分⽅程组。
微分⽅程中出现的未知函数的导数的最⾼阶解数称为微分⽅程的阶。
若⽅程中未知函数及其各阶导数都是⼀次的,称为线性常微分⽅程,⼀般表⽰为)()(')()(1)1(1)(t b y t a y t a y t a y n n n n =++++--若上式中的系数a i (t), i =1,2,…,n 均与t ⽆关,称之为常系数。
用MATLAB求解微分方程及微分方程组
dx dt X x , 0 即: Y y dy dt w dx ( X x) dt 2 2 ( X x ) (Y y ) 消去λ 得: dy w (Y y ) 2 2 dt ( X x ) (Y y )
1 0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8 -1 0 2 4 6 8 10 12
解
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);
用MATLAB求解微分方程
1. 微分方程的解析解
求微分方程(组)的解析解命令: dsolve(‘方程1’, ‘方程2’,…‘方程n’, ‘初始条件’, ‘自变量’)
记号: 在表达微分方程时,用字母 D 表示求微分,D2、D3 等 表示求高阶微分.任何 D 后所跟的字母为因变量,自变量可以指 定或由系统规则选定为确省. 例如,微分方程
任取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
(2)
(3)
3.因乙舰以速度v0沿直线x=1运动,设v0=1,则w=5,X=1,Y=t
matlab求解微分方程解析解
matlab求解微分方程解析解在数学和工程学科中,微分方程是一种重要的数学工具,它涉及到很多实际问题的模型和解决方法。
而Matlab作为一款强大的数学软件,可以方便地求解微分方程的解析解。
Matlab中求解微分方程的一种常见方法是使用符号计算工具箱(Symbolic Math Toolbox),它可以处理符号表达式和符号函数,包括微积分、代数、矩阵、符号等数学操作。
首先,我们需要定义微分方程的符号变量和初值条件。
例如,我们假设要求解的微分方程为dy/dx = x^2,初值条件为y(0)=1,则可以使用如下代码:syms x yode = diff(y,x) == x^2;cond = y(0) == 1;然后,我们可以将微分方程和初值条件作为参数传递给dsolve函数来求解微分方程的解析解。
例如:sol = dsolve(ode, cond);其中,sol为求解得到的符号表达式,可以使用vpa函数将其转换为数值解。
例如:sol_num = vpa(sol, 5);这样,我们就得到了微分方程的解析解,并将其转换为5位有效数字的数值解。
除了使用符号计算工具箱,Matlab还提供了许多数值方法来求解微分方程的数值解。
例如,可以使用ode45函数来求解微分方程的数值解。
例如,求解dy/dx = x^2,y(0)=1的数值解可以使用如下代码:fun = @(x,y) x^2;[t,y] = ode45(fun, [0,1], 1);其中,fun为微分方程的函数句柄,[0,1]为求解区间,1为初值条件。
t和y分别为求解得到的时间序列和解向量。
总之,Matlab提供了多种方法来求解微分方程的解析解和数值解,可以根据实际问题的需要选择不同的方法来求解。
实验二MATLAB数值计算常微分方程的求解
实验二MATLAB数值计算常微分方程的求解引言:微分方程是描述物理问题和工程问题的重要工具,也是数学和工程学科中的重要课题。
解析方法可以求得一些简单微分方程的解析解,但对于复杂问题,常常无法得到解析解。
此时,数值求解方法成为一种有效的工具。
MATLAB是一个功能强大的数值计算软件,对于求解常微分方程组也提供了多种数值方法。
本实验将介绍MATLAB中求解常微分方程(组)的方法,并通过实例来演示这些方法的应用。
一、MATLAB的常微分方程(组)求解函数MATLAB提供了多种函数用于求解常微分方程(组),其中最常用的函数是ode45、ode23和ode15s。
这些函数采用不同的数值方法,精度和效率也不同。
下面分别对这些函数进行简单介绍:1. ode45函数:ode45函数采用了一种自适应的步长控制算法,可以在一个时间段内自动调整步长。
这个函数的语法是:[t,y] = ode45(odefun,tspan,y0,options)其中,odefun是一个函数句柄,表示待解的常微分方程组,tspan是求解区间,y0是初始条件,options是一个结构体,包含其他的参数和选项。
2. ode23函数:ode23函数也采用了自适应的步长控制算法,但相对于ode45,它是一个简化版本,只允许使用比较简单的问题。
这个函数的语法与ode45相似:[t,y] = ode23(odefun,tspan,y0,options)3. ode15s函数:ode15s函数采用了一种隐式数值方法,适用于比较刚性(stiff)的常微分方程。
这个函数的语法与前两个函数也相似:[t,y] = ode15s(odefun,tspan,y0,options)以上是MATLAB提供的三种解常微分方程(组)的函数,根据问题的特点和求解的需求选择相应的函数。
二、实例演示在本实验中,我们将通过一个实例来演示如何使用MATLAB解常微分方程组。
假设有一个简单的线性常微分方程组:dy1/dt = -2y1 + y2dy2/dt = y1 - 2y2给定初始条件为y1(0)=1,y2(0)=0,求在t=[0,5]时,y1和y2的解。
实验二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为工具,探讨如何求解微分方程并对比解析解与数值解的差异。
一、引言微分方程是描述自然界中许多现象和问题的数学语言,它包含了未知函数及其导数与自变量之间的关系。
微分方程的求解可以帮助我们了解问题的性质和变化规律,并为实际应用提供参考。
在许多情况下,微分方程的解析解很难求得,这时我们可以利用计算机进行数值求解。
二、微分方程的数值解法1.欧拉法欧拉法是最简单的数值求解微分方程的方法之一。
它通过将微分方程转化为差分方程,然后利用离散的点逼近连续的解。
具体步骤如下:(1)将微分方程转化为差分方程,即用近似的导数代替真实的导数;(2)选择初始条件,即确定初始点的值;(3)选择步长和求解区间,即确定求解的范围和步长;(4)使用迭代公式计算下一个点的值;(5)重复步骤(4),直到达到指定的求解区间。
2.改进的欧拉法欧拉法存在精度较低的问题,为了提高精度,可以使用改进的欧拉法。
改进的欧拉法是通过使用两次导数的平均值来计算下一个点的值,从而提高了数值解的精度。
3.龙格-库塔法龙格-库塔法是一种常用的数值求解微分方程的方法,它通过使用多个点的导数来逼近连续解。
龙格-库塔法的步骤如下:(1)选择初始条件和步长;(2)使用迭代公式计算下一个点的值;(3)计算下一个点的导数;(4)根据导数的值和步长计算下一个点的值;(5)重复步骤(3)和(4),直到达到指定的求解区间。
龙格-库塔法的精度较高,适用于求解一阶和高阶微分方程。
三、微分方程的解析解解析解是指能够用公式或函数表示的方程的解。
有些微分方程具有解析解,可以通过数学方法求得。
例如,一阶线性常微分方程和某些特殊类型的二阶微分方程等。
解析解的优势在于精确性和直观性,能够帮助我们深入理解问题的本质。
如何使用MATLAB求解微分方程学习资料
如何使用MATLAB求解微分方程学习资料MATLAB是一种强大的数值计算和科学编程平台,可以用于求解微分方程和微分方程组。
在使用MATLAB求解微分方程之前,需要掌握一些基础知识,包括MATLAB的基本语法和常用的求解微分方程的技术。
下面是一些学习资料和步骤,帮助您使用MATLAB求解微分方程。
1.学习MATLAB基本语法和操作:首先,您需要学习MATLAB的基本语法和常用操作。
您可以参考MATLAB的官方文档、教程和手册,以及MATLAB的在线资源和视频教程。
这些资源可以帮助您掌握MATLAB的基本操作,建立良好的编程习惯。
2.学习求解微分方程的方法:在使用MATLAB求解微分方程之前,您需要了解一些常用的求解微分方程的方法,例如数值方法和解析方法。
数值方法包括欧拉法、龙格-库塔法、四阶龙格-库塔法等;解析方法包括分离变量法、线性微分方程的常系数齐次法和非齐次法等。
您可以参考微积分的教科书、在线资源和视频教程,掌握这些方法。
3. 使用MATLAB求解一阶微分方程:一阶微分方程是最简单的微分方程形式。
您可以首先尝试使用MATLAB求解一阶微分方程。
MATLAB提供了几个函数来求解一阶微分方程,例如ode45、ode23、ode113等。
您可以使用这些函数来解决特定的一阶微分方程,并观察结果。
可以使用plot 函数绘制微分方程的解,以获得更直观的理解。
4.使用MATLAB求解高阶微分方程:一旦您熟悉了使用MATLAB求解一阶微分方程的方法,您可以尝试使用同样的方法来求解高阶微分方程。
在求解高阶微分方程时,您需要将其转化为一组一阶微分方程。
例如,对于二阶线性微分方程,您可以引入一个新的变量来表示未知函数的导数,然后将其转化为一组一阶微分方程。
然后,您可以使用相同的求解函数来求解这组一阶微分方程。
5. 使用MATLAB求解微分方程组:对于多元微分方程组,MATLAB提供了更多的函数来求解。
例如,ode45s可以用于求解刚体动力学方程,ode23t可以用于求解刚体动力学方程。
matlab求微分方程精确解与近似解
Euler 折线法
考虑一维经典初值问题
dy dx
f (x, y) ,
y( x0 )
y0
,
x [a, b]
基本思想:用差商代替微商
根据 Talyor 公式,y(x) 在点 xk 处有
y( x) y( xk ) ( x xk ) y '( xk ) O(x2 )
y( xk1 ) y( xk ) hy '( xk ) O(h2 ) h xk1 xk
clear; global mu;mu=7; y0=[1;0]; [t,x]=ode45('verderpol',[0,40],y0); plot(t,x(:,1),'r-');
Matlab 求解微分方程小结
Matlab 函数
求解析解(通解或特解),用 dsolve 求数值解(特解),用 ode45、ode23 ...
Matlab 编程
Euler 折线法 Runga-Kutta 方法
上机作业
教材 P97:练习 1、2、3、4、5、6、7 要求:
请在上机之前将程序写好(参考附录),上机时 直接输入或修改附录程序即可。
例
4:求初值问题
dy dx
2 y 2x2
2x
的数值解,求解范
围为 [0,0.5]
y(0) 1
fun=inline('-2*y+2*x^2+2*x','x','y'); [x,y]=ode23(fun,[0,0.5],1);
注:也可以在 tspan 中指定对求解区间的分割,如:
[x,y]=ode23(fun,[0:0.1:0.5],1); %此时 x=[0:0.1:0.5]
matlab求微分方程组的解析解
MATLAB求微分方程组的解析解引言在科学与工程领域,微分方程组是一种常见的数学模型,用于描述各种物理现象和工程问题。
解析解是指能够用公式表达出来的精确解。
本文将介绍如何使用M ATL A B求解微分方程组的解析解。
1.方程组的建立首先,我们需要确定待求解的微分方程组。
假设我们有一个由n个微分方程组成的方程组,可以写为如下形式:d y1/dt=f1(t,y1,y2,...,yn)d y2/dt=f2(t,y1,y2,...,yn)......d y n/dt=f n(t,y1,y2,...,yn)其中`t`是自变量,`y1,y2,...,y n`是因变量,`f1,f2,...,fn`是给定的函数关系。
我们的目标是求解`y1(t),y2(t),...,yn(t)`的解析解。
2.使用MAT LAB求解M A TL AB提供了强大的求解微分方程组的工具,我们可以使用其中的函数来求解上述方程组的解析解。
首先,我们需要在MA T LA B中定义方程组的函数形式。
可以通过定义一个函数或者使用匿名函数来实现。
例如,我们可以定义一个名为`m yE qu at io ns`的函数,其输入参数为`t`和一个向量`y`,输出为一个向量`d y`,代表方程组的左侧和右侧的变量分别。
函数示例如下:f u nc ti on dy=m yE qua t io ns(t,y)%定义方程组d y=z er os(n,1);d y(1)=f1(t,y(1),y(2),...,y(n));d y(2)=f2(t,y(1),y(2),...,y(n));......d y(n)=fn(t,y(1),y(2),...,y(n));e n d然后,我们可以使用M AT LA B的`d so lv e`函数来求解微分方程组的解析解。
示例如下:s y ms ty1(t)y2(t)...yn(t)a s su me(t,'re al')a s su me(y1(t),'rea l')a s su me(y2(t),'rea l')......a s su me(y n(t),'rea l')e q n1=d if f(y1(t),t)==f1(t,y1(t),y2(t),...,y n(t));e q n2=d if f(y2(t),t)==f2(t,y1(t),y2(t),...,y n(t));......e q nn=d if f(yn(t),t)==fn(t,y1(t),y2(t),...,y n(t));e q ns=[eq n1,e qn2,...,eq nn];S=ds ol ve(e qn s);`S`即为方程组的解析解集合。
MATLAB解微分方程
0
500
1000 1500 2000 2500 3000
数学建模实验项目2(9)
五、实验题目
实验要求:针对下列问题,先建立数学模型,再确定模型 中的参数,并有结论。
1、求下列微分方程的解析解
(1) y′ − 2 y = x2
(2) y′ − xy = x, y(0) = 1
(3) y′′ + y′ = x, y(0) = 1, y′(0) = 1
x2
)
dx dt
−
x
=
0
x(0) = 2; x'(0) = 0
解: 令 y1=x,y2=y1’
则微分方程变为一阶微分方程组:
y2
'
=
y1'= 1000(1−
y2 y12
)y2
−
y1
y1(0) = 2, y2(0) = 0
(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);
其中k1是药物量被吸收到血液中的速率系数,k是血液中向体外 排除的速率系数,D是刚开始胃中或肌肉中的药物总量。 试用欧拉公式求上述微分方程数值解,并画出图形。(设
k1=0.6,k=0.2,D=200)
用Matlab软件求解微分方程 的解析解和数值解
数学与信息科学学院 孔祥庆
数学建模实验项目2(1)
一、实验名称:用Matlab软件求解微分方程的解析解和数值解 二、实验目的:
掌握用Matlab软件求解微分方程模型的解析解和数值解的方法 三、实验内容
matlab解微积分方程
matlab解微积分方程使用Matlab解微积分方程微积分方程是数学中的重要概念,广泛应用于物理学、工程学、经济学等领域。
解微积分方程是研究微积分方程的一个重要问题,而Matlab作为一种强大的数值计算软件,可以有效地解决微积分方程。
Matlab提供了多种求解微积分方程的方法,包括欧拉法、龙格-库塔法、四阶龙格-库塔法等。
这些方法可以用来求解常微分方程、偏微分方程以及一些特殊类型的微积分方程。
我们来看看如何使用Matlab求解常微分方程。
常微分方程是一种只涉及一个自变量的微分方程,可以表示为dy/dx = f(x, y),其中f(x, y)是已知的函数。
在Matlab中,可以使用ode45函数来求解常微分方程。
下面以一个简单的一阶常微分方程为例,来演示如何使用Matlab求解。
假设我们要求解方程dy/dx = x + y,且初始条件为y(0) = 1。
首先,我们需要定义方程的函数形式,即f(x, y) = x + y。
然后,使用ode45函数来求解:```function dydx = myode(x, y)dydx = x + y;end[t, y] = ode45(@myode, [0, 1], 1);```上述代码中,myode函数定义了方程的函数形式,ode45函数用于求解微分方程,[0, 1]表示求解的时间范围,1表示初始条件。
最后,得到的结果存储在变量t和y中,t表示时间,y表示方程的解。
除了常微分方程,Matlab还可以求解偏微分方程。
偏微分方程是一种涉及多个自变量的微分方程,可以表示为∂u/∂t = f(x, y, t, u, ∂u/∂x, ∂u/∂y)。
在Matlab中,可以使用pdepe函数来求解偏微分方程。
假设我们要求解一个简单的二维热传导方程,即∂u/∂t = ∂^2u/∂x^2 + ∂^2u/∂y^2,且初始条件为u(x, y, 0) = sin(x)sin(y),边界条件为u(0, y, t) = 0,u(π, y, t) = 0,u(x, 0, t) = 0,u(x, π, t) = 0。
matlab 求解微分方程
matlab 求解微分方程摘要:1.Matlab 简介2.微分方程基本概念3.Matlab 求解微分方程的方法4.常见微分方程求解实例5.总结正文:一、Matlab 简介Matlab 是一种广泛应用于科学计算、数据分析和可视化的编程语言。
它具有丰富的函数库和强大的矩阵计算能力,使得用户可以方便地完成各种复杂的数学运算和分析任务。
在微分方程求解领域,Matlab 同样具有很高的应用价值。
二、微分方程基本概念微分方程是数学中的一个重要分支,它描述了自然界和社会现象中许多变化规律。
微分方程可以分为偏微分方程和常微分方程两大类。
求解微分方程是数学和工程领域中的一个重要课题,关乎许多实际问题的解决。
三、Matlab 求解微分方程的方法Matlab 求解微分方程主要依赖于其内置的符号计算函数和数值计算函数。
用户可以根据微分方程的性质选择适当的求解方法,如符号解法、数值解法等。
Matlab 提供了丰富的函数和工具箱来支持微分方程的求解,如ode45、ode23 等。
四、常见微分方程求解实例1.常微分方程:例如一阶常微分方程y" + p(x)y = q(x),Matlab 可以通过ode45 函数求解。
2.偏微分方程:例如二维热传导方程,Matlab 可以通过pdepeye 函数求解。
3.线性微分方程组:例如常系数线性微分方程组,Matlab 可以通过ode45 等函数求解。
4.非线性微分方程:例如Riccati 方程,Matlab 可以通过ode45 等函数求解。
五、总结Matlab 作为一种强大的科学计算工具,可以帮助用户方便地求解各种微分方程。
matlab算法-求解微分方程数值解和解析解
MATLAB是一种用于数学计算、工程和科学应用程序开发的高级技术计算语言和交互式环境。
它被广泛应用于各种领域,尤其在工程和科学领域中被用于解决复杂的数学问题。
微分方程是许多工程和科学问题的基本数学描述,求解微分方程的数值解和解析解是MATLAB算法的一个重要应用。
1. 求解微分方程数值解在MATLAB中,可以使用各种数值方法来求解微分方程的数值解。
其中,常见的方法包括欧拉法、改进的欧拉法、四阶龙格-库塔法等。
这些数值方法可以通过编写MATLAB脚本来实现,从而得到微分方程的近似数值解。
以常微分方程为例,可以使用ode45函数来求解微分方程的数值解。
该函数是MATLAB中用于求解常微分方程初值问题的快速、鲁棒的数值方法,可以有效地得到微分方程的数值解。
2. 求解微分方程解析解除了求解微分方程的数值解外,MATLAB还可以用于求解微分方程的解析解。
对于一些特定类型的微分方程,可以使用符号计算工具箱中的函数来求解微分方程的解析解。
通过符号计算工具箱,可以对微分方程进行符号化处理,从而得到微分方程的解析解。
这对于研究微分方程的性质和特点非常有帮助,也有助于理论分析和验证数值解的准确性。
3. MATLAB算法应用举例在实际工程和科学应用中,MATLAB算法求解微分方程问题非常常见。
在控制系统设计中,经常需要对系统的动态特性进行分析和设计,这通常涉及到微分方程的建模和求解。
通过MATLAB算法,可以对系统的微分方程进行数值求解,从而得到系统的响应曲线和动态特性。
另外,在物理学、生物学、经济学等领域的建模和仿真中,也经常需要用到MATLAB算法来求解微分方程问题。
4. MATLAB算法优势相比于其他数学软件和编程语言,MATLAB在求解微分方程问题上具有明显的优势。
MATLAB提供了丰富的数值方法和工具,能够方便地对各种微分方程进行数值求解。
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 子目录,也可以保存在自定义文件夹,这时注意要增加搜索路径(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+uu。
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、符号说明:①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);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 +=+=plot(t1,p1,'-',t2,p2,'*')0510150.10.20.30.40.50.60.70.8图中‘*’曲线为战争中鲨鱼所占比例。