1. 程式人生 > >MATLAB求解常微分方程:ode45函式與dsolve函式

MATLAB求解常微分方程:ode45函式與dsolve函式

ode45函式無法求出解析解,dsolve可以求出解析解(若有),但是速度較慢.

1.      ode45函式

①求一階常微分方程的初值問題

[t,y] = ode45(@(t,y)y-2*t/y,[0,4],1);

plot(t,y);

求解 y’ – y + 2*t / y且初值y(0) = 1的常微分方程初值問題,返回自變數和函式的若干個值.

若不寫返回值,則會自動作出函式隨自變數的變化影象.

ode45(@(t,y)y-2*t/y,[0,4],1);

②求解一階微分方程組

x’ = -x^3-y,x(0)=1

y’ = x-y^3,y(0)=0.5.

自變數為t,且0<t<30.

求解過程如下.

第一步,在M函式檔案中將函式x和函式y寫成向量形式.

function f = fun(t,x);

f(1) = -x(1)^3 – x(2);

f(2) = x(1) – x(2)^3;

f = f(:);%確保f為列向量.

第二步,在M指令碼檔案中求解.

[t,x] = ode45(@fun,[0,30],[1;0.5]);

subplot(1,2,1);plot(t,x(:,1),t,x(:,2),':');xlabel('t');ylabel('x/y');%作x和y隨t變化圖

subplot(1,2,2);plot(x(:,1),x(:,2));xlabel('x');ylabel('y');%作x和y的相點陣圖

第三步,在命令視窗執行M指令碼檔案中的程式碼.

③求解高階常微分方程組

將高階常微分方程組通過變數替換轉化為一階的常微分方程組,然後用ode45求解.

2.      dsolve函式

①求解析解

y’ = a*x + b;

s = dsolve('D2y=a*y+b*x','x');

D2y用以表示y的二階導數,預設是以t為自變數的,所以最好指明自變數為x.

②初值問題

y’ = y – 2*t / y , y(0) = 1;

s = dsolve('Dy == y - 2*t / y','y(0) ==1');

③邊值問題

x*y’’ – 3*y’ = x^2 , y(1) = 0 , y(5) = 0;

s = dsolve('x*D2y - 3*Dy ==x^2','y(1)=0','y(5) == 0','x');

函式最後一個引數指明自變數為x.

④高階方程

求解y’’ = cos(2x) – y , y(0) = 1 , y’(0) = 0;

s=dsolve('D2y == cos(2*x) - y','y(0) =1','Dy(0) = 0','x');

simplify(s);

⑤方程組問題

f’ = f + g , g’ = -f + g,f(0) = 1, g(0) =2;

[f,g]= dsolve('Df == f + g','Dg = -f + g','f(0)==1','g(0) == 2','x');