Matlab在常微分方程求解中应用课件

Matlab在常微分方程求解中的应用Matlab在常微分方程求解中的应用1实验目的

（1）学会用Matlab软件求解微分方程的初值问题（2）了解微分方程数值解思想，掌握基本的微分方程数值解方法（3）学会根据实际问题建立简单微分方程数学模型（4）了解计算机数据仿真、数据模拟的基本方法

实验目的217世纪：初等解法18世纪：初等解法和无穷级数方法19世纪：解的存在性、奇点理论、定性理论、稳定性理论

包含一个自变量和它的未知函数以及未知函数的导数的等式

形成和发展与力学、天文学、物理学及其他自然科学技术的发展互相促进和推动

常微分方程17世纪：初等解法包含一个自变量和它的未知函数以及未知函数3定理

设函数在区域上连续，且在区域D内满足李普希兹(Lipschitz)条件，即存在正数L，使得对于R内任意两点与,恒有则初值问题(1)的解存在并且唯一。定理设函数在区域上连续，且在区域D内满足李普希兹(Li4常微分方程的解析解

求微分方程（组）的解析解命令:dsolve(‘方程1’,‘方程2’,…‘方程n’,‘初始条件’,‘自变量’)

结果：u=tg(t-c)记号：在表达微分方程时，用字母D表示求微分，D2、D3等表示求高阶微分。D后所跟字母为因变量，自变量可以指定或由系统规则选定为缺省。例如：微分方程可以表示为D2y=0.常微分方程的解析解求微分方程（组）的解析解命令:dsolv5

解输入命令:y=dsolve('D2y+4*Dy+29*y=0','y(0)=0,Dy(0)=15','x')结果为:y=3e-2xsin（5x）解输入命令:y=dsolve('D2y+4*Dy+26解输入命令：

[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=(c1-c2+c3+c2e-3t-c3e-3t)e2t

y=(c1e-4t+c2e-4t+c2e-3t-c3e-3t+c1-c2+c3)e2tz=(-c1e-4t+c2e-4t+c1-c2+c3)e2t

解输入命令：结果为：x=(c1-c2+c3+c27虽然说解析解是最精确的，但是实际问题中常要求研究常微分方程的数值解虽然说解析解是最精确的，但是实际问题中常要求研究常微分方程的8常微分方程的数值解

常微分方程中只有一些典型方程能求出初等解(用初等函数表示的解)。另外，有些初值问题虽然有初等解，但由于形式太复杂不便于应用。因此，有必要探讨常微分方程初值问题的数值解法。以下主要介绍一阶常微分方程初值问题几种经典数值解方法：欧拉法及改进的欧拉法。其它方法：龙格-库塔法、阿达姆斯方法；一阶微分方程组与高阶方程初值问题的数值解法；二阶常微分方程值问题的差分方法等。常微分方程的数值解常微分方程中只有一些典型方程9求解常微分方程初值问题的数值解的整体思路：

（1）寻求准确解在一系列离散节点：上的近似值

称为问题的数值解，数值解所满足的离散方程统称为差称为步长，实用中常取定步长。分格式，求解常微分方程初值问题的数值解的整体思路：（1）寻求准10建立数值解法的一些途径1、用差商代替导数

若步长h较小，则有故有公式：此即欧拉法建立数值解法的一些途径1、用差商代替导数若步长h较小，则11Euler法的几何意义:找到了积分的一条近似折线Euler法的几何意义:找到了积分的一条近似折线122、使用数值积分对方程y`=f(x,y),两边由xi到xi+1积分，并利用梯形公式，有：实际应用时，与欧拉公式结合使用：此即改进的欧拉法故有公式：2、使用数值积分对方程y`=f(x,y),两边由xi到xi13Matlab在常微分方程求解中应用课件14Matlab在常微分方程求解中应用课件15Matlab在常微分方程求解中应用课件163、使用泰勒公式

以此方法为基础，有龙格-库塔法、线性多步法等方法。4、数值公式的精度

当一个数值公式的截断误差可表示为O（hk+1）时（k为正整数，h为步长），称它是一个k阶公式。k越大，则数值公式的精度越高。欧拉法是一阶公式，改进的欧拉法是二阶公式。龙格-库塔法有二阶公式和四阶公式。线性多步法有四阶阿达姆斯外插公式和内插公式。3、使用泰勒公式以此方法为基础，有龙格-库塔法、线性17ODE指令列表

MATLAB用於求解常微分方程式的指令:

指令方法应用ODE类别ode45ExplicitRunge-Kutta(4,5)pairofDormand-PrinceNonstiffODEode23ExplicitRunge-Kutta(2,3)pairofBogackiandShampineNonstiffODEode113VariableorderAdams-Bashforth-MoultonPECEsolverNonstiffODEode15sNumericaldifferentiationformulas(NDFS)StiffODEode23sModifiedRosenbrockformulaoforder2StiffODEode23tTrapezoidalrulewitha“free”interpolantStiffODEode23tbImplicitRunge-KuttaformulawithabackwarddifferentiationformulaofordertwoStiffODEODE指令列表MATLAB用於求解常微分方程式的指令:18

适用于

Nonstiff系統

速率（即微分值）差异相常大使用一般的ode45、ode23或ode113來求解,可能会使得积分的步长(StepSizes)变得很小,以便降低积分误差至可容忍范围以內,会导致计算时间过长专门对付Stiff系统的指令，例如ode15s、ode23s、

ode23t及ode23tb指令项目繁多，最主要可分两大类：适用于

Nonstiff系統

一般的常微分方程式都是Nonstiff系統直接采用ode45、ode23或ode113來求解指令项目繁多，最主要可分两大类：适用于Nonstiff系19提示使用Simulink來求解常微分方程式Simulink是和MATLAB共同使用的一套软件可使用拖拉的方式來建立动力系统可直接产生C语言代码或进行动画演示功能非常强大提示使用Simulink來求解常微分方程式20ODE指令基本用法

使用ODE指令时，必须先将要求解的ODE表示

成一个函数

输入为t(时间)及y(状态函数:StateVariables）输出则为dy（状态变量的微分值）ODE指令基本用法使用ODE指令时，必须先将要求解的OD21ODE指令基本用法

若ODE函数的文件为odeFile.m,则调用ODE指令

的格式如下：

[t,y]=solver('odeFile',[t0,t1],y0)

[t0,t1]是积分的时间区间y0代表起始条件（InitialConditions）solver是前表所列的各种ODE指令t是输出的时间向量y是对应的状态变量向量ODE指令基本用法若ODE函数的文件为odeFile.m22[t,y]=solver（’f’,ts,y0,options）ode45ode23ode113ode15sode23s由待解方程写成的m-文件名ts=[t0,tf]，t0、tf为自变量的初值和终值函数的初值ode23：组合的2/3阶龙格-库塔-芬尔格算法ode45：运用组合的4/5阶龙格-库塔-芬尔格算法自变量值函数值用于设定误差限(缺省时设定相对误差10-3,绝对误差10-6),命令为：options=odeset（’reltol’,rt,’abstol’,at）,rt，at：分别为设定的相对误差和绝对误差.[t,y]=solver（’f’,ts,y0,options23

1、在解n个未知函数的方程组时，y0和y均为n维向量，m-文件中的待解方程组应以y的分量形式写成.

2、使用Matlab软件求数值解时，高阶微分方程必须等价地变换成一阶微分方程组.注意:2、使用Matlab软件求数值解时，高阶微分方程必须等价24VanderPol微分方程1928年荷兰的范德波耳(VanderPol)为描述LC回路的电子管振荡器建立了著名的vanderPol方程.它在自激振荡理论中有着重要的意义,一直作为数学物理方程中的一个基本方程.这是一个具有可变非线性阻尼的微分方程,代表了一类极为典型的非线性问题.和其他非线性微分方程在数学上无法精确求解一样,人们一直在努力寻找求解这类方程近似解析解的方法,并乐于用VanderPol方程来检验求解方法的有效性.

VanderPol微分方程1928年荷兰的范德波耳(Va25VanderPol微分方程其方程形式为：

令

将VanderPol微分方程化成标准形式VanderPol微分方程其方程形式为：将Vander26写成向量的形式：

为一个向量，代表状态变量写成向量的形式：为一个向量，代表状态变量27设

=1，ODE文件(vdp1.m)显示如下：

>>typevdp1.m

functiondy=vdp1(t,y)mu=1;dy=[y(2);mu*(1-y(1)^2)*y(2)-y(1)];建立了vdp1.m后，即可选用前述ODE指令来求解设=1，ODE文件(vdp1.m)显示如下：28在=1时，vanderPol方程并非Stiff系统，所以使用ode45来显示结果

odeBasic01.mode45('vdp1',[025],[33]');

[025]代表积分的时间区间，[33]’则代表起始条件因为没有输出变量，所以上述程序执行结束后，MATLAB只画出状态变量对时间的图形在=1时，vanderPol方程并非Stiff系统29Matlab在常微分方程求解中应用课件30

odeGetData01.m[t,y]=ode45('vdp1',[025],[33]');plot(t,y(:,1),t,y(:,2),':');xlabel('Timet');ylabel('Solutiony(t)andy''(t)');legend('y(t)','y''(t)');[t,y]=ode45('vdp1',[025],31Matlab在常微分方程求解中应用课件32也可以画出及在相位平面（PhasePlane）的运动情况

odePhastPlot01.m[t,y]=ode45('vdp1',[025],[33]');plot(y(:,1),y(:,2),'-o');xlabel('y(t)');ylabel('y''(t)');也可以画出及在相位平面（PhasePlane）33Matlab在常微分方程求解中应用课件34当值越来越大时，VanderPol方程就渐变为一个Stiff系统，此时若要求解该类问题，就必须用对应于Stiff系统的指令

当值越来越大时，VanderPol方程就渐变为一个S35将值改成1000，ODE文件改写成（vdp2.m）：

>>typevdp2.mfunctiondy=vdp2(t,y)mu=1000;dy=[y(2);mu*(1-y(1)^2)*y(2)-y(1)];将值改成1000，ODE文件改写成（vdp2.36选用专门用来解决Stiff系统问题的ODE指令，例如ode15s，來求解此系统并作图显示ode15s01.m[t,y]=ode15s('vdp2',[03000],[21]');subplot(2,1,1);plot(t,y(:,1),'-o');xlabel('Timet');ylabel('y(t)');subplot(2,1,2);plot(t,y(:,2),'-o');xlabel('Timet');ylabel('y''(t)'); 选用专门用来解决Stiff系统问题的ODE指令，例如37由上图可知，的变化相当剧烈（超过），这也就是Stiff系统的特点Matlab在常微分方程求解中应用课件38绘制二维平面相位图：ode15s02.m若要产生在某些特定时间点的状态变量值，则调用ODE指令的格式可改成：

[t,y]=solver('odeFile',[t0,t1,…,tn],y0)其中[t0,t1,…,tn]即是特定时间点所形成的向量[t,y]=ode15s('vdp2',[03000],[21]');subplot(1,1,1);plot(y(:,1),y(:,2),'-o');xlabel('y(t)');ylabel('y''(t)')绘制二维平面相位图：[t,y]=ode15s('vdp239Matlab在常微分方程求解中应用课件40ODE指令的选项ODE指令可以接受第四个输入变量，代表积分过程用到的各种选项（Options），此种ODE指令的格式为：

[t,y]=solver('odeFile',[t0,tn],y0,options)其中options是由odeset指令來控制的结构变量结构变量即包含了积分过程用到的各种选项odeset的一般格式如下：options=odeset('name1',value1,'name2',value2,…)name1的值为value1，name2的值为value2…

ODE指令的选项ODE指令可以接受第四个输入变量，代表积41也可以只改变一个现有的options中某栏位的值，格式如下：newOptions=odeset(options,'name',value);若要获取某栏位的值，可用odeget，格式如下：

value=odeget(otpions,'name');当使用odeset指令時，若无输入任何变量名，则odeset会显示所有变量名及其值，并以大括号代表预设值也可以只改变一个现有的options中某栏位的值，格式如下：42>>odesetAbsTol:[positivescalarorvector{1e-6}]RelTol:[positivescalar{1e-3}]NormControl:[on|{off}]NonNegative:[vectorofintegers]OutputFcn:[function_handle]OutputSel:[vectorofintegers]Refine:[positiveinteger]Stats:[on|{off}]InitialStep:[positivescalar]MaxStep:[positivescalar]BDF:[on|{off}]MaxOrder:[1|2|3|4|{5}]Jacobian:[matrix|function_handle]JPattern:[sparsematrix]Vectorized:[on|{off}]Mass:[matrix|function_handle]MStateDependence:[none|{weak}|strong]MvPattern:[sparsematrix]MassSingular:[yes|no|{maybe}]InitialSlope:[vector]Events:[function_handle]>>odeset43由odeset产生的ODE选项

类别栏位名称资料形态预设值说明誤差容忍度之相关栏位RelTol正数量相对误差容忍度AbsTo1正数量或向量绝对误差容忍度积分输出之相关栏位OutPutFcn字串‘odeplot’输出函式（若ODE指令无输出变量，則在数值积分执行完毕后，MATLAB会呼叫此输出函数）OutputSel索引向量全部ODE指令之输出变量的索引值，以決定那些输出变量之元素将被送到输出函式Refine正整数1或4(forode45)Refine=2可产生两倍数量的输出点，Refine=3可产生三倍数量的输出点，依此类推Statson或offoffStats=‘on’产生计算过程的各种统计资料由odeset产生的ODE选项类别栏位名称资料形态44常用到的档位说明

f在积分误差容忍度方面，每一次积分所产生的局部误差e(i)，必须满足下列方程式：

max(RelTol*,AbsTol(i))其中i代表第i个状态

降低RelTol及AbsTol来求得更精确的积分结果常用到的档位说明f在积分误差容忍度方面，每一次积分所产生的45odeRelTol01.msubplot(2,1,1);ode45('vdp1',[025],[33]');title('RelTol=0.01');options=odeset('RelTol',0.00001);subplot(2,1,2);ode45('vdp1',[025],[33]',options);title('RelTol=0.0001');odeRelTol01.msubplot(2,1,1);46第一个图所使用的相对误差值是0.01（预设值），第二个图所使用的相对误差值是0.00001，因此我们得到较细密的点，但所花的计算时间也会比较长

第一个图所使用的相对误差值是0.01（预设值），第二个图所使47Matlab在常微分方程求解中应用课件48积分输出说明积分输出的相关处理

选用一个OutputFcnOutputFcn的预设值是“odeplot”，其功能为画出所有的状态变量的图像其它可用的函式odephas2：画出2-D的相位平面（PhasePlane）odephas3：画出3-D的相位平面odeprint：随时在指令视窗显示计算结果积分输出说明积分输出的相关处理49options=odeset('OutputFcn','odephas2');%使用odephas2进行绘图ode45('vdp1',[025],[33]',options);options=odeset('OutputFcn',50Matlab在常微分方程求解中的应用Matlab在常微分方程求解中的应用51实验目的

（1）学会用Matlab软件求解微分方程的初值问题（2）了解微分方程数值解思想，掌握基本的微分方程数值解方法（3）学会根据实际问题建立简单微分方程数学模型（4）了解计算机数据仿真、数据模拟的基本方法

实验目的5217世纪：初等解法18世纪：初等解法和无穷级数方法19世纪：解的存在性、奇点理论、定性理论、稳定性理论

包含一个自变量和它的未知函数以及未知函数的导数的等式

形成和发展与力学、天文学、物理学及其他自然科学技术的发展互相促进和推动

常微分方程17世纪：初等解法包含一个自变量和它的未知函数以及未知函数53定理

设函数在区域上连续，且在区域D内满足李普希兹(Lipschitz)条件，即存在正数L，使得对于R内任意两点与,恒有则初值问题(1)的解存在并且唯一。定理设函数在区域上连续，且在区域D内满足李普希兹(Li54常微分方程的解析解

求微分方程（组）的解析解命令:dsolve(‘方程1’,‘方程2’,…‘方程n’,‘初始条件’,‘自变量’)

结果：u=tg(t-c)记号：在表达微分方程时，用字母D表示求微分，D2、D3等表示求高阶微分。D后所跟字母为因变量，自变量可以指定或由系统规则选定为缺省。例如：微分方程可以表示为D2y=0.常微分方程的解析解求微分方程（组）的解析解命令:dsolv55

解输入命令:y=dsolve('D2y+4*Dy+29*y=0','y(0)=0,Dy(0)=15','x')结果为:y=3e-2xsin（5x）解输入命令:y=dsolve('D2y+4*Dy+256解输入命令：

[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=(c1-c2+c3+c2e-3t-c3e-3t)e2t

y=(c1e-4t+c2e-4t+c2e-3t-c3e-3t+c1-c2+c3)e2tz=(-c1e-4t+c2e-4t+c1-c2+c3)e2t

解输入命令：结果为：x=(c1-c2+c3+c257虽然说解析解是最精确的，但是实际问题中常要求研究常微分方程的数值解虽然说解析解是最精确的，但是实际问题中常要求研究常微分方程的58常微分方程的数值解

常微分方程中只有一些典型方程能求出初等解(用初等函数表示的解)。另外，有些初值问题虽然有初等解，但由于形式太复杂不便于应用。因此，有必要探讨常微分方程初值问题的数值解法。以下主要介绍一阶常微分方程初值问题几种经典数值解方法：欧拉法及改进的欧拉法。其它方法：龙格-库塔法、阿达姆斯方法；一阶微分方程组与高阶方程初值问题的数值解法；二阶常微分方程值问题的差分方法等。常微分方程的数值解常微分方程中只有一些典型方程59求解常微分方程初值问题的数值解的整体思路：

（1）寻求准确解在一系列离散节点：上的近似值

称为问题的数值解，数值解所满足的离散方程统称为差称为步长，实用中常取定步长。分格式，求解常微分方程初值问题的数值解的整体思路：（1）寻求准60建立数值解法的一些途径1、用差商代替导数

若步长h较小，则有故有公式：此即欧拉法建立数值解法的一些途径1、用差商代替导数若步长h较小，则61Euler法的几何意义:找到了积分的一条近似折线Euler法的几何意义:找到了积分的一条近似折线622、使用数值积分对方程y`=f(x,y),两边由xi到xi+1积分，并利用梯形公式，有：实际应用时，与欧拉公式结合使用：此即改进的欧拉法故有公式：2、使用数值积分对方程y`=f(x,y),两边由xi到xi63Matlab在常微分方程求解中应用课件64Matlab在常微分方程求解中应用课件65Matlab在常微分方程求解中应用课件663、使用泰勒公式

以此方法为基础，有龙格-库塔法、线性多步法等方法。4、数值公式的精度

当一个数值公式的截断误差可表示为O（hk+1）时（k为正整数，h为步长），称它是一个k阶公式。k越大，则数值公式的精度越高。欧拉法是一阶公式，改进的欧拉法是二阶公式。龙格-库塔法有二阶公式和四阶公式。线性多步法有四阶阿达姆斯外插公式和内插公式。3、使用泰勒公式以此方法为基础，有龙格-库塔法、线性67ODE指令列表

MATLAB用於求解常微分方程式的指令:

指令方法应用ODE类别ode45ExplicitRunge-Kutta(4,5)pairofDormand-PrinceNonstiffODEode23ExplicitRunge-Kutta(2,3)pairofBogackiandShampineNonstiffODEode113VariableorderAdams-Bashforth-MoultonPECEsolverNonstiffODEode15sNumericaldifferentiationformulas(NDFS)StiffODEode23sModifiedRosenbrockformulaoforder2StiffODEode23tTrapezoidalrulewitha“free”interpolantStiffODEode23tbImplicitRunge-KuttaformulawithabackwarddifferentiationformulaofordertwoStiffODEODE指令列表MATLAB用於求解常微分方程式的指令:68

适用于

Nonstiff系統

速率（即微分值）差异相常大使用一般的ode45、ode23或ode113來求解,可能会使得积分的步长(StepSizes)变得很小,以便降低积分误差至可容忍范围以內,会导致计算时间过长专门对付Stiff系统的指令，例如ode15s、ode23s、

ode23t及ode23tb指令项目繁多，最主要可分两大类：适用于

Nonstiff系統

一般的常微分方程式都是Nonstiff系統直接采用ode45、ode23或ode113來求解指令项目繁多，最主要可分两大类：适用于Nonstiff系69提示使用Simulink來求解常微分方程式Simulink是和MATLAB共同使用的一套软件可使用拖拉的方式來建立动力系统可直接产生C语言代码或进行动画演示功能非常强大提示使用Simulink來求解常微分方程式70ODE指令基本用法

使用ODE指令时，必须先将要求解的ODE表示

成一个函数

输入为t(时间)及y(状态函数:StateVariables）输出则为dy（状态变量的微分值）ODE指令基本用法使用ODE指令时，必须先将要求解的OD71ODE指令基本用法

若ODE函数的文件为odeFile.m,则调用ODE指令

的格式如下：

[t,y]=solver('odeFile',[t0,t1],y0)

[t0,t1]是积分的时间区间y0代表起始条件（InitialConditions）solver是前表所列的各种ODE指令t是输出的时间向量y是对应的状态变量向量ODE指令基本用法若ODE函数的文件为odeFile.m72[t,y]=solver（’f’,ts,y0,options）ode45ode23ode113ode15sode23s由待解方程写成的m-文件名ts=[t0,tf]，t0、tf为自变量的初值和终值函数的初值ode23：组合的2/3阶龙格-库塔-芬尔格算法ode45：运用组合的4/5阶龙格-库塔-芬尔格算法自变量值函数值用于设定误差限(缺省时设定相对误差10-3,绝对误差10-6),命令为：options=odeset（’reltol’,rt,’abstol’,at）,rt，at：分别为设定的相对误差和绝对误差.[t,y]=solver（’f’,ts,y0,options73

1、在解n个未知函数的方程组时，y0和y均为n维向量，m-文件中的待解方程组应以y的分量形式写成.

2、使用Matlab软件求数值解时，高阶微分方程必须等价地变换成一阶微分方程组.注意:2、使用Matlab软件求数值解时，高阶微分方程必须等价74VanderPol微分方程1928年荷兰的范德波耳(VanderPol)为描述LC回路的电子管振荡器建立了著名的vanderPol方程.它在自激振荡理论中有着重要的意义,一直作为数学物理方程中的一个基本方程.这是一个具有可变非线性阻尼的微分方程,代表了一类极为典型的非线性问题.和其他非线性微分方程在数学上无法精确求解一样,人们一直在努力寻找求解这类方程近似解析解的方法,并乐于用VanderPol方程来检验求解方法的有效性.

VanderPol微分方程1928年荷兰的范德波耳(Va75VanderPol微分方程其方程形式为：

令

将VanderPol微分方程化成标准形式VanderPol微分方程其方程形式为：将Vander76写成向量的形式：

为一个向量，代表状态变量写成向量的形式：为一个向量，代表状态变量77设

=1，ODE文件(vdp1.m)显示如下：

>>typevdp1.m

functiondy=vdp1(t,y)mu=1;dy=[y(2);mu*(1-y(1)^2)*y(2)-y(1)];建立了vdp1.m后，即可选用前述ODE指令来求解设=1，ODE文件(vdp1.m)显示如下：78在=1时，vanderPol方程并非Stiff系统，所以使用ode45来显示结果

odeBasic01.mode45('vdp1',[025],[33]');

[025]代表积分的时间区间，[33]’则代表起始条件因为没有输出变量，所以上述程序执行结束后，MATLAB只画出状态变量对时间的图形在=1时，vanderPol方程并非Stiff系统79Matlab在常微分方程求解中应用课件80

odeGetData01.m[t,y]=ode45('vdp1',[025],[33]');plot(t,y(:,1),t,y(:,2),':');xlabel('Timet');ylabel('Solutiony(t)andy''(t)');legend('y(t)','y''(t)');[t,y]=ode45('vdp1',[025],81Matlab在常微分方程求解中应用课件82也可以画出及在相位平面（PhasePlane）的运动情况

odePhastPlot01.m[t,y]=ode45('vdp1',[025],[33]');plot(y(:,1),y(:,2),'-o');xlabel('y(t)');ylabel('y''(t)');也可以画出及在相位平面（PhasePlane）83Matlab在常微分方程求解中应用课件84当值越来越大时，VanderPol方程就渐变为一个Stiff系统，此时若要求解该类问题，就必须用对应于Stiff系统的指令

当值越来越大时，VanderPol方程就渐变为一个S85将值改成1000，ODE文件改写成（vdp2.m）：

>>typevdp2.mfunctiondy=vdp2(t,y)mu=1000;dy=[y(2);mu*(1-y(1)^2)*y(2)-y(1)];将值改成1000，ODE文件改写成（vdp2.86选用专门用来解决Stiff系统问题的ODE指令，例如ode15s，來求解此系统并作图显示ode15s01.m[t,y]=ode15s('vdp2',[03000],[21]');subplot(2,1,1);plot(t,y(:,1),'-o');xlabel('Timet');ylabel('y(t)');subplot(2,1,2);plot(t,y(:,2),'-o');xlabel('Timet');ylabel('y''(t)'); 选用专门用来解决Stiff系统问题的ODE指令，例如87由上图可知，的变化相当剧烈（超过），这也就是Stiff系统的特点Matlab在常微分方程求解中应用课件88绘制二维平面相位图：ode15s02.m若要产生在某些特定时间点的状态变量值，则调用ODE指令的格式可改成：

[t,y]=solver('odeFile',[t0,t1,…,tn],y0)其中[t0,t1,…,tn]即是特定时间点所形成的向量[t,y]=ode15s('vdp2',[03000],[21]');subplot(1,1,1);plot(y(:,1),y(:,2),'-o');xlabel('y(t)');ylabel('y''(t)')绘制二维平面相位图：[t,y]=ode15s('vdp289Matlab在常微分方程求解中应用课件90ODE指令的选项ODE指令可以接受第四个输入变量，代表积分过程用到的各种选项（Options），此种ODE指令的格式为：

[t,y]=solver('odeFile',[t0,tn],y0,options)其中options是由odeset指令來控制的结构变量结构变量即包含了积分过程用到的各种选项odeset的一般格式如下：options=odeset('name1',value1,'name2',value2,…)name1的值为value1，name2的值为value2…

ODE指令的选项ODE指令可以接受第四个输入变量，代表积91也可以只改变一个现有的options中某栏位的值，格式如下：newOptions=odeset(options,'name',value);若要获取某栏位的值，可用odeget，格式如下：

value=odeget(otpions,'name');当使用odeset指令時，若无输入任何变量名，则odeset会显示所有变量名及其值，并以大括号代表预设值也可以只改变一个现有的options中某栏位的值，格式如下：92>>odesetAbsTol:[positivescalarorvector{1e-6}]RelTol:[positivescalar{1e-3}]NormControl:[on|{off}]NonNegative:[vectorofintegers]OutputFcn:[function_handle]OutputSel:[vectorofintegers]Refine:[positiveinteger]Stats:[on|{off}]InitialStep:[pos