用MATLAB求解非线性微分方程

总结一下MATLAB中求解微分方程的思路和步骤。固然，网上很多关于此类的技术型文章，但往往一看下来发现，文章中的友情链接比文章字数还多，要了解这一篇文章，你要先了解那个；要了解那个，你又要了解那个那个；以此类推。了解完了，花都谢了。我在这里要写的是娱乐型文章，即试图用最少的字数陈述出MATLAB的求解微分方程方法，看完后可以记得清楚且不亦乐乎。其实个人十分推荐MATLAB中的帮助文件，内容实在是特别给力，简明的不能再简明了。可惜了是英文写的，即使是天天要面对英语的我，也颇为头疼。那么以下我将结合我对MATLAB中帮助文档的理解，提出我对求解最简单微分方程的方法的总结。以求解一个经典的非线性方程为例。

         做一个最基本的假设：你们都看过高数。

 

一。老湿发话了：童鞋们，求解一下这个方程，判断她是否稳定。要是稳定，那么她是否存在极限环：

一看明白了，这不就是传说中的范德普方程。地球人都知道她稳定并有极限环。现在我们就看看如何用MATLAB求解她的轨迹。

 

二。一般的计算机求解方程的方法无外乎是这样：首先把该方程改写成一个规范的形式，一般使用状态空间表示法；而后调用已有的算法进行求解；最后对得出的结果进行处理，比如画图之类的。接下来就对这三大步分别作出解释。

 

三。输入待求解的方程。

         首先我们知道，状态空间的标准形式(自由系统)是：

这里X是列向量，F是作用于列向量的函数，可以是线性也可以是非线性。范德普方程可以改写成这样的标准形式:

MATLAB中关于输入输入待解方程的语句特别简单。需要先定义一个普通函数，函数名任意，姑且叫做myFcn，格式如下

function xdot = myFcn (t, x)  

 

需要注意的是，函数必须含有t, x两个参数，名称可以自己任意定。xdot是这个函数本身的返回值，只出现在这个函数内部，因此也可以任意定。

定义中的x被视为是一个列向量，x(i)表示列向量中的第i个分量。那么F函数的每一个分量即简单地用表达时给出即可。其中的自变量可以引用x(i)。以范德普方程为例：

 

xdot = [x(2) ; u(1-x(1)^2)*x(2)-x(1)]  

 

于是，这两句话便构成了待解函数。

 

 

四。调用MATLAB函数进行求解

         通常人工求解微分方程需要知道初始值，计算机求解也不例外。另外，由于非线性方程一般只有数值解，故计算精度也可以调整。这些都是可以自己调整的参数。

         调用MATLAB计算求解常微分方程的模式很简单，格式为：

 

[t, x] = ode45(@myFcn, [开始时间 结束时间], [初始值列向量], options)

 

注意到求解出来的结果是[t, x]即一堆数，所以需要我们进行后处理比如画图之类的。

ode45表示了MATLAB中的一种内置计算微分方程的算法，最为常用，出于娱乐目的，就先用这个。

@表示待求解的方程，如myFcn。

[开始时间 结束时间]表示求解的时间段，不必定义间隔。

[初始值列向量]就不用多说了。通常在求解之前定义一个变量x0列向量表示初始值，然后输入起来就方便多了。

options本身是一个变量。注意，她的名称不固定，但是他是一个以结构体为类型的变量。options很重要，稍后讨论。

 

五。进行后处理。在得到[t, x]后可以自己用plot神马的进行画图。不过我一般不这样，各位看官不必惊慌。

 

六。options的用法。

         options在MATLAB帮助中是这样定义格式的：

         options = odeset (“Name”, Value, “Name”, Value, …)

         意思是，options这个变量要用到odeset这个函数来对他进行赋值。而odeset这个函数的参数必须是成对出现的：一个名称对应一个数值。

         我们要用到的是这样的：

 

         options= odeset(“OutputFcn”, @odephas2, “reltol”, 1e-5);

 

注意，odeset中的参数名称不可任意定，因为都是预定义好的。”OutputFcn”使用预定义的函数进行输出，而与定义的函数有好几种，使用@进行选择，我们选odephas2即画相平面图(phase portrait)。之后的那个参数表示相对误差的最大允许值，不说也明白。

         有一个问题，用odephas2画出的图图不好看，因为上面打了好多圈圈。解决办法是找到该文件，修改其中所有的plot语句即可，把那个”-o”去掉就行了。

 

七。就是这个样了。总结一下，求解范德普方程可以用如下语句：

 

首先是函数定义：

function xdot = myFcn (t, x)

xdot = [x(2) ; u(1-x(1)^2)*x(2)-x(1)]

 

其次是主程序：

% Solving options predetermined

options=odeset('OutputFcn',@odephas2,'reltol',1e-5);

% Solving vdp equation withdifferent ini conditions

for a = -3 : 0.1 : 3   

   b1=sqrt(25-a^2);

   b2=-sqrt(25-a^2);

   % Solve the problem

   [t,y]=ode45(@myFcn, [0 20], [a b1], options);

   hold on

   [t,y]=ode45(@myFcn, [0 20] ,[a b2], options);

   hold on    

end

grid on