命令集

 [time,x]=solver(str,t,x0)   計算ODE或由字串str 給定的ODE的值，部分解已在向量time中給出。在向量time

中給出部分解，包含的是時間值。還有部分解在矩陣x中給出，x的列向量每

個方程在這些值下的解。對於標量問題，方程的解將在向量x中給出。這些解

在時間區間t(1)到t(2)上計算得到。其初始值是x0 即x(t(1)).此方程組有str 指定的

M檔案中函式表示出。這個函式需要兩個引數：標量t 和向量x, 應該返回向量

x’(即x的導數)。因為對標量ODE來說，x和x’都是標量。在M檔案中輸入odefile

可得到更多資訊。同時可以用命令numjac來計算Jacobi函式。

[t,x]=solver(str,t,x0,val)   此方程的求解過程同上，結構val包含使用者給solver的命令。參見odeset和表1，可得到更多資訊。

Ode45     此方法被推薦為首選方法。

Ode23     這是一個比ode45低階的方法。

Ode113     用於更高階或大的標量計算。

Ode23t     用於解決難度適中的問題。

Ode23s     用於解決難度較大的微分方程組。對於系統中存在常量矩陣的情況也有用。

Ode15s     與ode23相同，但要求的精度更高。

Ode23tb     用於解決難度較大的問題，對於系統中存在常量矩陣的情況也有用。

Set=odeset(set1,vak1,set2,val2,…)  返回結構set,其中包含用於ODE求解方程的設定引數，   有關可用設定

的資訊參見表1。

Odeget(set,’set1’)    返回結構set中設定set1的值。

有許多設定對odeset控制的ODE解是有用的，參見表1。例如，如果要在求解過程中畫出解的圖形，可以

輸入：inst=odeset(‘outputfcn’,’odeplot’);.也可使用命令odedemo。

表1    ODE求解方程的設定引數

RelTol   給出求解方程允許的相對誤差

AbsTol  給出求解方程允許的絕對誤差

Refine  給出與輸入點數相乘的因子

OutputFcn 這是一個帶有輸入函式名的字串，該字串將在求解函式執行的每步被呼叫：odephas2(畫出2D

的平面相點陣圖)。Odephas3( 畫出3D 的平面相點陣圖)，odeplot( 畫出解的圖形)，odeprint( 顯示中間結

果)

OutputSel 是一個整型向量。指出哪些元素應該被傳遞給函式，特別是傳遞給OutputFcn

Stats  如果引數Stats為on,則將統計並顯示出計算過程中資源消耗情況

Jacobian… 如果編寫ODE檔案程式碼以便F（t,y,jocobian）返回dF/dy,則將jacobian設定為on

Jconstant…如果雅可比數df/dy是常量，則將此引數設定為on

Jpattern…  如果編寫ODE檔案的編碼以便函式F（[],[],jpattern）返回帶有零的稀疏矩陣並輸出非零元素dF/fy,則

需將Jpattern設定為on

Vectorized…如果編寫ODE檔案的編碼以便函式F(t,[y1,y2……])返回[F(t,y1)  F(t,y2)…],則將此引數設定成on

Events…   如果ODE檔案中帶有引數‘events’，則將此引數設定成on

Mass…   如果編寫ODE檔案編碼以實現函式F(t,[],‘mass’)返回M和M(t),應將此引數設定成on

MassConstant…如果矩陣M(t)是常量，則將此引數設定成on

MaxStep…   此引數是限定演算法能使用的區間長度上限的標量

InitialStep…   給出初始步長的標量。如果給定的區間太大，演算法就使用一個較小的步長

MaxOrder…   此引數只能被ode15s使用，它主要是指定ode15s的最高階數,並且此引數應是從1到5的整數

BDF…   此引數只能被ode15s使用，如果倒推微分公式而不是用通常所使用的微分公式，則要將它設定為

on

NormControl…如果演算法根據norm(e)<=max(Reltol*norm(y),Abstol)來步積分過程中的錯誤，則要將它設定成on

例1

求解x’=-x^2 初值x(0)=1

建立函式xprim1,將此函式儲存在M檔案xprim1.m中：

function xprim=xprim1(t,x)

xprim=-x.^2;

然後呼叫MATLAB的ODE演算法求解方程。然後畫出解的圖形：

[t,x]=ode45(‘xprim1’,[0,1],1);

plot(t,x,’-‘,t,x,’o’);

xlabel(‘time t0=0,tt=1’);

ylabel(‘x values x(0)=1’);

或者

首先建立xprim2,將此函式儲存在M檔案xprim2.m中：

function xprim=xprim2(t,x)

xprim=x.^2;

然後呼叫MATLAB的ODE演算法求解方程。然後畫出解的圖形：

[t,x]=ode45(‘xprim2’,[0,0.95],1);

注意：在MATLAB中計算出的點在微分絕對值大的區域內更密集些。

如果想在影象中直接觀察，則在程式碼後輸入以下命令

plot(t,x,’o‘,t,x,’-’);

xlabel(‘time t0=0,tt=0.95’);

ylabel(‘x values x(0)=1’);

例2求解下列二元一階微分方程組

X1’=X1-0.1X1X2+0.01t

X2’=-X2+0.02X1X2+0.044

X1(0)=30

X2(0)=20

這個方程組用在人口動力學中。可以認為是單一化的捕食者---被捕食者模式。例如，狐狸和兔子。表示被捕食者， 表示捕食者。如果被捕食者有無限的食物，並且不會出現捕食者。於是有，這個式子是以指數形式增長的。大量的被捕食者將會使捕食者的數量增長；同樣，越來越少的捕食者會使被捕食者的數量增長。而且，人口數量也會增長。

建立xprim3,將此函式儲存在M檔案xprim3.m 中：

function xprim=xprim3(t,x)

xprim=[x(1)-0.1*x(1)*x(2)+0.01*t; …

-x(2)+0.02*x(1)*x(2)+0.04*t];

然後呼叫一個ODE演算法和畫出解的圖形：

[t,x]=ode45(‘xprim3’,[0,20],[30;20]);

plot(t,x);

xlabel(‘time t0=0,tt=20’);

ylabel(‘x values x1(0)=30,x2(0)=20’);

在MATLAB中,也可以根據x2函式繪製出x1圖形,命令plot(x(:2),x(:1)) 可繪製出平面相點陣圖.

例3 帶引數的微分方程

X1’=a-(b+1)X1+X1^2X2

X2’=bX1+X1^2X2

X1(0)=1

X2(0)=3

方程由下面的M檔案stiff1.m定義:

functionstiff=stiff1(t,x)

global a;   %變數不能放入引數表中

glabol b;

stiff=[0,0];    %stiff必須是一個冒號變數

stiff(1)=a-(b+1)*x1+x(1)^2*x(2);

stiff(2)=b*x(1)-x(1)^2*x(2);

下面的M檔案給出一個比較困難的問題：

global a;a=100;

global b;b=1;

tic;

[t,X]=ode23(‘stiff1’，[0 10],[1 3]);

toc

size(t)

執行後得到的結果如下：

elapsed_time=72.1647

ans=34009

用專門解決複雜問題的解法ode23s,將得到較好的結果：

elapsed_time=1.0098

ans=103

對於邊界值問題，除了微分方程，還有邊界處的值。在一維下這意味著至少有兩個條件。

例4假設要研究一根杆的溫度分佈情況。這根杆一端的溫度是T0,另一端的溫度是T1。令y(x)表示這根杆的溫度，函式f(x)表示加熱源。從時間t=0開始，在相當長的時間內加熱這根杆，直至達到平衡狀態。這就是所謂的定常值或穩定狀態。這個定常值可由下面的方程模型表示：

-Y’(X)=f(X),0<X<1

Y(0)=T0

Y(1)=T1

假設這根杆兩端為：x=0和x=1。

假設在其兩端又一根固定的柱子（或者可以看成是一個連線兩個島嶼的橋），如圖7所示。

令y(x)表示載入函式g(x)後彎曲的柱子。此問題需要有兩個關於此柱子兩端的邊界條件。假設這根柱子非常牢

固的固定在牆上，即y在牆上的導數是0。可以得到下面的ODE,其中介紹了自然協調系統：

Y’’’(X)=g(x),0<X<1

Y(0)= Y’(0)= Y(1)= Y’(1)=0

解y(x)的導數可由有限的差分代替

如果用這些差分方程來代替ODE中的導數，就能得到一個所有未知的yi的方程組。其係數矩陣是一個有序區

間，此區間的寬度決定於這個微分方程的導數個數.

根據前面的溫度模型的方程研究一下杆的溫度分佈，將所有的導數換成不同的差分並得到：

其中fj=f(xj)。為了簡單起見，設M=6,即給定的y0和y6，而y1,y2,….y5 為未知變數。於是就有

 

注意，y0=T0和yM=T1必須移到方程組的右邊。此時得到的矩陣是一個對角矩陣，其對角線上的元素為2，並且上一對角線和下一對角線上的元素為1。

下面解此問題的檔案temperature.m。使用者必須先給出分段數及f(x)(用點符號)，最後給出T0和T1。有關稀疏矩陣的更多資訊參見其他資料。

%杆上的溫度分佈，用T0和T1分別表示兩端溫度
%這根杆放在x座標的0和1 區間上，並被分成M個子區間，每個子區間的長度為1/M
%建立稀疏矩陣方程Ax=b並求解
%矩陣A對角陣，並以稀疏矩陣的形式儲存
clear;
M=input(‘Give the number ofsubintervals (M):’);
Deltax=1/M;
xx=0:deltax:1;
funcStr= input(‘give f(x),the extra heatsource(e.g.,x.^3):’,’s’);
T0=input(‘Give y(0) (left): ’);
T1=input(‘Give y(1) (right): ’);
%構造對角陣和方程右邊b
vectorOnes=ones(M-1,1);
A=spdiags([-vectorOnes,2*vectorOnes,-vectorOnes],[-1 0 1],M-1,M-1);
x=xx(2:end-1);  %x為區域內的值。
f=eval(funcStr);  %響應的f(x)的值。
b=deltax^2f;
b(1)=b(1)+T0; %對邊界值x=0,x=1 進行特殊處理。
b(end)=b(end)+T1;
b=b’;
%解線性方程
y=A\b;  %y在區間內：j=1:M-1.
y=[T0;y;T1];  %y在整個區間內：0<=x<=1.
clf;
%上面圖形表示外部熱源。
%下面圖形表示杆上的熱分佈。
subplot(2,1,1);
plot(x,f);
grid on;
title(‘External heatsource f(x).’,’FontSize’,14);
subplot(2,1,2);
plot(xx,y,’r’);
grid on;
title(‘Tempearture distribution in a rod.’,’FontSize’,14);
%將區間分成等份，根據方程f(x)=在圖中可以得到解。

例5 Matlab 中二階常微分方程的數值解法

y''-3y'+2y=1;  y(0)=1; y'(0)=0; 用數值解法求y(0.5)

[email protected](t,x)[x(2);3*x(2)-2*x(1)+1];

[t,y]=ode45(odefun,[0:0.01:2],[1 0]);

plot(t,y)

[t y]

結果

y(0.5000)=0.7896

y= dsolve('D2y-3*Dy+2*y=1','Dy(0)=0','y(0)=1');

>> y

y =exp(t) - exp(2*t)/2 + 1/2

>> feval(@(t)exp(t) - exp(2*t)/2 +1/2,0.5)

ans = 0.7896

例6求解範德普方程

首先是函式定義：
 
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

例7x1.^2-10x1+x2.^2+b=0

     x1*x2.^2+x1-10x2+8=0

用MATLAB實現，編寫非線性方程組的M檔案fx1.m如下所示：

   function y=fx1(x)

  y(1)=x(1)*x(1)-10*x(1)+x(2)*x(2)+8;

  y(2)=x(1)*x(2)*x(2)+x(1)-10*x(2)+8;

   y=;

編寫非線性方程組導數的M檔案dfx1.m

  function y=dfx1(x)

     y(1)=2*x(1)-10;

     y(2)=2*x(2);

     y(3)=x(2)*x(2)+1;

     y(4)=2*x(1)*x(2)-10;

     y=;

附一篇文章《用MATLAB求解非線性微分方程》

原始地址：http://blog.csdn.net/anhuixue/article/details/7560558

總結一下MATLAB中求解微分方程的思路和步驟。MATLAB中的幫助檔案，內容實在是特別給力，簡明的不能再簡明瞭。可惜了是英文寫的，即使是天天要面對英語的我，也頗為頭疼。那麼以下我將結合我對MATLAB中幫助文件的理解，提出我對求解最簡單微分方程的方法的總結。以求解一個經典的非線性方程為例。

一求解一下這個方程，判斷她是否穩定。要是穩定，那麼她是否存在極限環;

二.一般的計算機求解方程的方法：

首先把該方程改寫成一個規範的形式，一般使用狀態空間表示法；

呼叫已有的演算法進行求解；

最後對得出的結果進行處理，比如畫圖之類的。

三.輸入待求解的方程。

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”去掉就行了。