[轉] http://blog.sina.com.cn/s/blog_46e9b2010100tsqv.html
用matlab時間也不短了，可是一直沒有接觸過微分方程。這次看看書，學習學習，記點兒筆記。
| y=dsolve(f1,f2,...,fmO; y=dsolve(f1,f2,...,fm,'x'); |
|
syms t;
u=exp(-5*t)*cos(2*t-1)+5; uu=5*diff(u,t,2)+4*diff(u,t)+2*u; syms t y; y=dsolve(['D4y+10*D3y+35*D2y+50*Dy+24*y=87*exp(-5*t)*cos(2*t-1)+92*exp(-5*t)*sin(2*t-1)+10'])
yc=latex(y)
|
將yc的內容copy到latex中編譯，得到結果。
關於Matlab的微分方程，直到今天才更新第2篇，實在是很慚愧的事——因為原因都在於太懶惰，而不是其他的什么。
在上一篇中，我們使用dsolve可以解決一部分能夠解析求解的微分方程、微分方程組，但是對於大多數微分方程（組）而言不能得到解析解，這時數值求解也就是沒有辦法的辦法了，好在數值解也有很多的用處。
數值分析方法中講解了一些Eular法、 Runge-Kutta 法等一些方法，在matlab中內置的ode求解器可以實現不同求解方法的相同格式的調用，而不必太關心matlab究竟是用什么算法完成的。
這一回我們來說明ode45求解器的使用方法。
1.ode45求解的上手例子：
求解方程組
Dx=y+x(1-x^2-y^2);
Dy=-x+y*(1-x^2-y^2)
初值x=0.1;y=0.2;
先說明一下最常用的ode45調用方式，和相應的函數文件定義格式。
[t,x]=ode45(odefun,tspan,x0);
其中，Fun就是導函數，tspan為求解的時間區間（或時間序列，如果采用時間序列，則必須單調），x0為初值。
這時，函數文件可以采用如下方式定義
function dx=odefun(t,x)
對於上面的小例子，可以用如下的程序求解。
| function jixianhuan
function dx=jxhdot(t,x) |
2.終值問題
tspan可以是遞增序列，也可以為遞減序列，若為遞減則可求解終值問題。
[t,x]=ode45(@zhongzhiode,[3,0],[1;0;2]);plot(t,x)
function dx=zhongzhiode(t,x)
dx=[2*x(2)^2-2;
-x(1)+2*x(2)*x(3)-1;
-2*x(2)+2*x(3)^2-4];
結果如下
3.odeset
options = odeset('name1',value1,'name2',value2,...)
[t,x]=solver(@fun,tspan,x0,options)
通過odeset設置options
第一，通過求解選項的設置可以改善求解精度，使得原本可能不收斂的問題收斂。
options=odeset('RelTol',1e-10);
第二，求解形如M(t,x)x'=f(t,x)的方程。
例如，方程
x'=-0.2x+yz+0.3xy
y'=2xy-5yz-2y^2
x+y+z-2=0
可以變形為
[1 0 0][x'] [-0.2x+yz+0.3xy]
[0 1 0][y'] = [2xy-5yz-2y^2 ]
[0 0 0][z'] [x+y+z-2 ]
這樣就可以用如下的代碼求解該方程
function mydae
M=[1 0 0;0 1 0;0 0 0];
options=odeset('Mass',M);
x0=[1.6,0.3,0.1];
[t,x]=ode15s(@daedot,[0,1.5],x0,options);plot(t,x)
function dx=daedot(t,x)
dx=[
-0.2*x(1)+x(2)*x(3)+0.3*x(1)*x(2);
2*x(1)*x(2)-5*x(2)*x(3)-2*x(2)*x(2);
x(1)+x(2)+x(3)-2];
4.帶附加參數的ode45
有時我們需要研究微分方程組中的參數對於解的影響，這時采用帶有參數的ode45求解會使求解、配合循環使用，可以使得求解的過程更加簡捷。
使用方法：只需將附加參數放在options的后面就可以傳遞給odefun了。
看下面的例子。
function Rossler
clear;clc
a=[0.2,0.2];
b=[0.2,0.5];
c=[5.7,10];
x0=[0 0 0];
for jj=1:2
[t,x]=ode45(@myRossler,[0,100],x0,[],a(jj),b(jj),c(jj));
figure;plot3(x(:,1),x(:,2),x(:,3));grid on;
end
function dx=myRossler(t,x,a,b,c)
dx=[
-x(2)-x(3);
x(1)+a*x(2);
b+(x(1)-c)*x(3)];
5. 剛性方程的求解
剛性方程就是指各個自變量的變化率差異很大，會造成通常的求解方法失效。
這是matlab中自帶的一個例子，使用ode15s求解，如果用ode45求解就會出現錯誤。
function myode15study
[t,Y] = ode15s(@vdp1000,[0 3000],[2 0]);
plot(T,Y(:,1),'-o')
figure;plot(Y(:,1),Y(:,2))
function dy = vdp1000(t,y)
dy = zeros(2,1);
dy(1) = y(2);
dy(2) = 1000*(1 - y(1)^2)*y(2) - y(1);
6.高階微分方程的求解
通常的方法是進行變量替換，將原方程降階，轉換成更多變量的一階方程組進行求解。
在這個例子里我們求解一個動力學系統里最常見的一個運動方程
,其中f=sin(t)
function myhighoder
clear;clc
x0=zeros(6,1);
[t,x]=ode45(@myhigh,[0,100],x0);
plot(t,x(:,1))
function dx=myhigh(t,x)
f=[sin(t);0;0];;
M=eye(3);
C=eye(3)*0.1;
K=eye(3)-0.5*diag(ones(2,1),1)-0.5*diag(ones(2,1),-1);
dx=[x(4:6);inv(M)*(f-C*x(4:6)-K*x(1:3))];
7.延遲微分方程
matlab提供了dde23求解非中性微分方程。dde23的調用格式如下：
sol = dde23(ddefun,lags,history,tspan)
lags是延遲量，比如方程中包含y1(t-0.2)和y2(t-0.3)則可以使用lags=[0.2,0.3]。
這里的ddefun必須采用如下的定義方式：
dydt = ddefun(t,y,Z)
其中的Z(:,1)就是y(t-lags(1)),Z(:,2)就是y(t-lags(2)）...
下面是個使用dde23求解延遲微分方程的例子。
function mydde23study
% The differential equations
%
% y'_1(t) = y_1(t-1)
% y'_2(t) = y_1(t-1)+y_2(t-0.2)
% y'_3(t) = y_2(t)
%
% are solved on [0, 5] with history y_1(t) = 1, y_2(t) = 1, y_3(t) = 1 for
% t <= 0.
clear;clc
lags=[1,0.2];
history=[1;1;1];
tspan=[0,5];
sol = dde23(@myddefun,lags,history,tspan)
plot(sol.x,sol.y)
function dy = myddefun(t,y,Z)
dy=[
Z(1,1);
Z(1)+Z(2,2);
y(2) ];
8.ode15i求解隱式微分方程
[T,Y] = ode15i(odefun,tspan,y0,yp0)
yp0為y'的初值。
odefun的格式如下 dy = odefun(t,y,yp),yp表示y'，而方程中應該使得f(t,y,y')=0
function myodeIMP
% The problem is
%
% y(1)' = -0.04*y(1) + 1e4*y(2)*y(3)
% y(2)' = 0.04*y(1) - 1e4*y(2)*y(3) - 3e7*y(2)^2
% y(3)' = 3e7*y(2)^2
%
% It is to be solved with initial conditions y(1) = 1, y(2) = 0, y(3) = 0
% to steady state.
clear;clc
y0=[1;0;0];
fixed_y0=[1;1;1];
yp0=[0 0 0];
fixed_yp0=[];
[y0mod,yp0mod]=decic(@myodefunimp,0,y0,fixed_y0,yp0,fixed_yp0);
tspan=[0, logspace(-6,6)];
[t,y] = ode15i(@myodefunimp,tspan,y0mod,yp0mod);
y(:,2)=1e4*y(:,2);
semilogx(t,y)
function res=myodefunimp(t,y,yp)
res=[
-yp(1)-0.04*y(1)+1e4*y(2)*y(3);
-yp(2)+0.04*y(1)-1e4*y(2)*y(3)-3e7*y(2)^2;
-yp(3)+3e7*y(2)^2;
];
這次要接觸一個新的求解ode的方法，就是使用simulink的積分器求解。
1.還是做我們研究過的一個例子（在初識matlab微分方程(2)中采用的）。
Dx=y+x(1-x^2-y^2);
Dy=-x+y*(1-x^2-y^2)
初值x=0.1;y=0.2;
積分器中設置初始條件；f(u)中指定Dx,Dy的計算公式。
運行這個仿真，scope中可以看到兩個變量的時程如下：
在WorkSpace里可以得到tout和yout，執行plot(yout(:,1),yout(:,1))得到與ode45求解相似的結果如下
2.這部分解決一個使用ode求解器dde23沒法求解的一類延遲微分方程(中性微分方程)。
形如x'(t)=f(x'(t-t1),x(t),x(t-t2),x(t-t3))這類方程。dde23是無法求解的，但是可以借助simulink仿真求解。
看下面的這個例子。
x'(t)=A1*x(t-t1)+A2*x'(t-t2)+B*u(t)
t1=0.15;t2=0.5
A1=[-12 3 -3] A2=[0.02 0 0] B=[0]
[106 -116 62] [0 0.03 0] [1]
[207 -207 113] [0 0 0.04] [2]
在continuous里找到transport Delay,就可以實現對於信號的延遲，因此可以建立如下仿真模型
從而在scope中可以得到如下仿真結果
OK~初識微分方程到了這里我想應該可以做個終結，因為我想作為零基礎的材料來看，到這里也就可以了。以后還可能再有微分方程的內容，還請感興趣的朋友多捧場吧。
最后，大力推薦一本書薛定宇老師的《高等應用數學問題的Matlab求解》，確實很經典。學習Matlab的時間也不算短了，可是每次翻看這本書總是能讓我有溫故而知新的感覺，是我目前見過的最好的Matlab書。強烈推薦!（對於從來沒有接觸過matlab的人來說或許有點兒難，但是如果你以后要用matlab的話買一本絕對不會后悔的。）