關(guān)于MATLAB連續(xù)模型求解方法介紹和分析

連續(xù)模型是指模型是連續(xù)函數(shù)的一類模型總稱，具體建模方法主要是微分方程建模。微分方程建模是數(shù)學(xué)建模的重要方法，因?yàn)樵S多實(shí)際問(wèn)題的數(shù)學(xué)描述將導(dǎo)致求解微分方程的定解問(wèn)題。把形形色色的實(shí)際問(wèn)題化成微分方程的定解問(wèn)題，大體上可以按以下幾步：

1.根據(jù)實(shí)際要求確定要研究的量（自變量、未知函數(shù)、必要的參數(shù)等）并確定坐標(biāo)系。

2.找出這些量所滿足的基本規(guī)律（物理的、幾何的、化學(xué)的或生物學(xué)的等等）。

3.運(yùn)用這些規(guī)律列出方程和定解條件。

MATLAB 在微分模型建模過(guò)程中的主要作用是求解微分方程的解析解， 將微分方程轉(zhuǎn)化為一般的函數(shù)形式。 另外， 微分方程建模， 一定要做數(shù)值模擬， 即根據(jù)方程的表達(dá)形式， 給出變量間關(guān)系的圖形， 做數(shù)值模擬也需要用 MATLAB 來(lái)實(shí)現(xiàn)。

微分方程的形式多樣，微分方程的求解也是根據(jù)不同的形式采用不同的方法， 在建模比賽中， 常用的方法有三種：

1.用 dsolve 求解常見的微分方程解析解

2.用 ODE 家族的求解器求解數(shù)值解

3.使用專用的求解器求解

1.常規(guī)微分方程的求解

微分方程在 MATLAB 中固定的表達(dá)方式， 這些基本的表達(dá)方式如下表所示：

對(duì)于通常的微分方程， 一般需要先求解析解， 那么 dsolve 是首先考慮的求解器，因?yàn)閐solve 能夠求解解析解，其具體的用法如下：

[實(shí)例]求微分方程 xy'+y-e^x=0 在初始條件 y(1)=2e 下的特解，并畫出解函數(shù)的圖形．

求解本問(wèn)題的 Matlab 程序?yàn)椋?/p>

syms x y

y=dsolve('x*Dy+y-exp(x)=0','y(1)=2*exp(1)','x')

ezplot(y)

微分方程的特解為：y=1/x*exp(x)+1/x*exp (1) (Matlab格式)，即 y=(e+e^x)/x ，此函數(shù)的圖形如圖 1：

圖1y關(guān)于x的函數(shù)圖象

2.ODE 家族求解器

如果微分方程的解析形式求解不出來(lái)， 那么退而求其次的辦法是求解數(shù)值解，那么這個(gè)時(shí)候就需要用 ODE 家族的求解器求解微分方程的數(shù)值解啦。

因?yàn)闆]有一種算法可以有效地解決所有的 ODE 問(wèn)題，為此，MATLAB 提供了多種求解器，對(duì)于不同的 ODE 問(wèn)題，采用不同的 Solver。MATLAB 中常用的微分方程數(shù)值解的求解器及特點(diǎn)如下表所示。

要特別提醒的是：ode23、ode45 是極其常用的用來(lái)求解非剛性標(biāo)準(zhǔn)形式一階常微分方程(組)初值問(wèn)題解的 Matlab 的常用程序，其中：

ode23 采用龍格-庫(kù)塔2 階算法，用3 階公式作誤差估計(jì)來(lái)調(diào)節(jié)步長(zhǎng)，具有低等的精度．

ode45 則采用龍格-庫(kù)塔4 階算法，用5 階公式作誤差估計(jì)來(lái)調(diào)節(jié)步長(zhǎng)，具有中等的精度．

[實(shí)例]導(dǎo)彈追蹤問(wèn)題

設(shè)位于坐標(biāo)原點(diǎn)的甲艦向位于 x 軸上點(diǎn) A(1, 0) 處的乙艦發(fā)射導(dǎo)彈，導(dǎo)彈頭始終對(duì)準(zhǔn)乙艦。如果乙艦以最大的速度v0(是常數(shù))沿平行于 y 軸的直線行駛，導(dǎo)彈的速度是 5*v0，求導(dǎo)彈運(yùn)行的曲線方程，以及乙艦行駛多遠(yuǎn)時(shí)，導(dǎo)彈將它擊中？

記導(dǎo)彈的速度為 w，乙艦的速率恒為v0。設(shè)時(shí)刻 t乙艦的坐標(biāo)為 (X(t),Y(t))，導(dǎo)彈的坐標(biāo)為 (x(t),y(t))。當(dāng)零時(shí)刻，(X(0),Y(0))=(1,0)，(x(0),y(0))=(0,0)，建立微分方程模型：

因乙艦以速度v0沿直線 x=1 運(yùn)動(dòng)，設(shè)v0=1，w=5，X=1，Y=t，因此導(dǎo)彈運(yùn)動(dòng)軌跡的參數(shù)方程為：

MATLAB 求解數(shù)值解程序如下：

(1)定義方程的函數(shù)形式：

function dy=eq2(t,y)

dy=zeros(2,1);

dy(1)=5*(1-y(1))/sqrt((1-y(1))^2+(t-y(2))^2);

dy(2)=5*(t-y(2))/sqrt((1-y(1))^2+(t-y(2))^2);

(2)求解微分方程的數(shù)值解

t0=0,tf=0.21;

[t,y]=ode45('eq2',[t0 tf],[0 0]);

X=1;Y=00.21;plot(X,Y,'-')

plot(y(:,1),y(:,2),'*'),hold on

x=01; y=-5*(1-x).^(4/5)/8+5*(1-x).^(6/5)/12+5/24;

plot(x,y,'r')

3.專用求解器

對(duì)于復(fù)雜的微分方程模型的求解， 可以借助 MATLAB 偏微分方程工具箱中的專用求解器。以下將以一個(gè)實(shí)例來(lái)看看如何借助偏微分方程工具箱來(lái)實(shí)現(xiàn)一個(gè)微分方程的求解與數(shù)值仿真。

所研究的對(duì)象是一個(gè)二階波的方程：

這個(gè)時(shí)候要查看一下 MALTAB 中哪個(gè)函數(shù)能求解相類似的方程， solvepde 可以求解的方程形式為：

可以發(fā)現(xiàn)只要通過(guò)參數(shù)設(shè)定就可以將所要求解的方程轉(zhuǎn)化成這種標(biāo)準(zhǔn)形式。

具體求解步驟如下：

(1) 設(shè)置參數(shù)

c = 1;

a = 0;

f = 0;

m = 1;

(2) 定義波的空間位置

numberOfPDE = 1;

model = createpde(numberOfPDE);

geometryFromEdges(model,@squareg);

pdegplot(model,'EdgeLabels','on');

ylim([-1.1 1.1]);

axis equal

title'Geometry With Edge Labels Displayed';

xlabel x

ylabel y

(3) 定義微分方程模型的系數(shù)和邊界條件

specifyCoefficients(model,'m',m,'d',0,'c',c,'a',a,'f',f);

applyBoundaryCondition(model,'dirichlet','Edge',[2,4],'u',0);

applyBoundaryCondition(model,'neumann','Edge',([1 3]),'g',0);

(4)定義該問(wèn)題的有限元網(wǎng)格

generateMesh(model);

figure

pdemesh(model);

ylim([-1.1 1.1]);

axis equal

xlabel x

ylabel y

(5)定義初始條件

u0 = @(location) atan(cos(pi/2*location.x));

ut0 = @(location) 3*sin(pi*location.x).*exp(sin(pi/2*location.y));

setInitialConditions(model,u0,ut0);

(6)方程的求解

n = 31;% 求解次數(shù)

tlist = linspace(0,5,n);

model.SolverOptions.ReportStatistics ='on';

result = solvepde(model,tlist);

u = result.NodalSolution;

(7)模型的數(shù)值仿真

figure

umax = max(max(u));

umin = min(min(u));

for i = 1:n

pdeplot(model,'XYData',u(:,i),'ZData',u(:,i),'ZStyle','continuous',...

'Mesh','off','XYGrid','on','ColorBar','off');

axis([-1 1 -1 1 umin umax]);

caxis([umin umax]);

xlabel x

ylabel y

zlabel u

M(i) = getframe;

end