1、求解一階常微分方程
d y d t = a y 2 \cfrac{dy}{dt}=ay^2 dtdy?=ay2
clc,clear
syms y(t) a %定義符號變量
dsolve(a*y^2-diff(y)==0)
結(jié)果
ans =
0
-1/(C1 + a*t)
2、求解三階常微分方程
d 3 y d t 3 = b y \cfrac{d^3y}{dt^3}=by dt3d3y?=by
clc,clear
syms y(t) b %定義符號變量
dsolve(diff(y,3)-b*y==0)
結(jié)果
ans =
C3*exp(b^(1/3)*t) + C1*exp(-t*((3^(1/2)*b^(1/3)*1i)/2 + b^(1/3)/2)) + C2*exp(t*((3^(1/2)*b^(1/3)*1i)/2 - b^(1/3)/2))
3.求解常微分方程
x 2 + y + ( x ? 2 y ) y ′ = 0 x^2+y+(x-2y)y\prime=0 x2+y+(x?2y)y′=0
clc,clear
syms y(x) %定義符號變量
dsolve(x^2+y+(x-2*y)*diff(y)==0)
結(jié)果：
ans =
x/2 + ((4*x^3)/3 + x^2 + C1)^(1/2)/2
x/2 - ((4*x^3)/3 + x^2 + C1)^(1/2)/2
4. 求解二階常微分方程
d 2 y d t 2 = a y , 初 始 條 件 為 y ( 0 ) = 5 , y ′ ( 0 ) = 1 \cfrac{d^2y}{dt^2}=ay, \text初始條件為y(0)=5,y\prime(0)=1 dt2d2y?=ay,初始條件為y(0)=5,y′(0)=1
clc,clear
syms y(t) a %定義符號變量
dy = diff(y);
y = dsolve(diff(y,2)==a*y,y(0)==5,dy(0)==1);
y = simplify(y) %化簡結(jié)果
結(jié)果
y =
(exp(a^(1/2)*t)*(5*a^(1/2) + 1))/(2*a^(1/2)) + (exp(-a^(1/2)*t)*(5*a^(1/2) - 1))/(2*a^(1/2))
5.求解常微分方程組
{ d y d t = z + 5 t d z d t = ? 2 y \left\{\begin{aligned} \cfrac{dy}{dt}&=z+5t \\ \\ \cfrac{dz}{dt}&=-2y \\ \end{aligned}\right. ??????????????dtdy?dtdz??=z+5t=?2y?
clc,clear
syms y(t) z(t) %定義符號變量
dy = diff(y);
dz = diff(z);
[y1,z1] = dsolve(dy==z+5*t,dz==-2*y);
y1 = simplify(y1) %化簡結(jié)果
z1 = simplify(z1) %化簡結(jié)果
結(jié)果：
y1 =
(2^(1/2)*C2*cos(2^(1/2)*t))/2 + (2^(1/2)*C1*sin(2^(1/2)*t))/2 + 5/2
z1 =
C1*cos(2^(1/2)*t) - 5*t - C2*sin(2^(1/2)*t)
6.求解初值為下述常微分方程的解析解和精確解
clc,clear
yx = @(x,y) -2*y+2*x^2+2*x; %定義微分方程右端的匿名函數(shù)
[x,y] = ode45(yx,[0,0.5],1)
[x1,y1] = ode23(yx,[0,0.5],1)
plot(x,y,x1,y1)
7.應(yīng)用兩個模型分別計(jì)算并預(yù)測未來五年的工資走勢
根據(jù)以下提供的某省份年平均工資統(tǒng)計(jì)表（單位：元），通過馬爾薩斯模型和阻滯增長模型構(gòu)建微分方程，求解并繪制圖像，應(yīng)用兩個模型分別計(jì)算并預(yù)測未來五年的工資走勢.
x = [1978:1:2010]
y=[566,632,745,755,769,789,985,1110,1313,1428,1782,1920,2150,2292,2601,3149,4338,5145,5809,6271,6854,7656,8772,10007,11374,12567,14332,16614,19228,22844,26404,29638,32074]
7.1 馬爾薩斯模型
模型假設(shè)：
- 設(shè)x(t)表示t時刻的工資，且x(t)連續(xù)可微
- 工資的增長率r是常數(shù)
- 工資的變化是封閉的，工資的變化只和工作的年限有關(guān)，工作年限越長，工資越高。
模型建立和求解：
假設(shè)t時刻到
t
+
Δ
t
t+\Delta t
t+Δt時間內(nèi)工資的增量為：
{
d
x
d
t
=
r
x
,
x
(
0
)
=
x
0
,
\left\{\begin{aligned} \cfrac{dx}{dt}&=rx, \\ \\ x(0)&=x_0, \\ \end{aligned}\right.
??????????dtdx?x(0)?=rx,=x0?,?
其解為：
x
(
t
)
=
x
0
?
e
r
t
x(t)=x_0*e^{rt}
x(t)=x0??ert
7.2 阻滯增長模型
對增長率r進(jìn)行修正。如果在工資較少的時候，把r看成常數(shù)，工資增長到一定數(shù)量的時候，就應(yīng)該把r看作應(yīng)該隨著工資增加而減少的量，即將增長率r表示成工資x(t)的函數(shù)r(x),r(x)是x的減函數(shù)。
模型假設(shè)：
- 設(shè)r(x)為x的線性函數(shù)，r(x)=r-sx(工程師原則，首選線性)
- 最大工資是: x m x_m xm?，即當(dāng)x= x m x_m xm?時，增長率r( x m x_m xm?)=0.
模型求解：
由假設(shè)1、2, 可得
r
(
x
)
=
r
(
1
?
x
x
m
r(x)=r(1-\cfrac{x}{x_m}
r(x)=r(1?xm?x?), 則有
{
d
x
d
t
=
r
(
1
?
x
x
m
)
x
,
x
(
t
0
)
=
x
0
,
\left\{\begin{aligned} \cfrac{dx}{dt}&=r(1-\cfrac{x}{x_m})x, \\ \\ x(t_0)&=x_0, \\ \end{aligned}\right.
??????????dtdx?x(t0?)?=r(1?xm?x?)x,=x0?,?
其解為：
x
(
t
)
=
x
m
1
+
(
x
m
x
0
?
1
)
e
?
r
(
t
?
t
0
)
x(t)=\cfrac{x_m}{1+(\cfrac{x_m}{x_0}-1)e^{-r(t-t_0)}}
x(t)=1+(x0?xm???1)e?r(t?t0?)xm??
x=[566 632 745 755 769 789 985 1110 ...
1313 1428 1782 1920 2150 2292 2601 3149 ...
4338 5145 5809 6271 6854 7656 8772 10007 ...
11374 12567 14332 16614 19228 22844 26404 29638 ...
32074]';
t=1978:2010;t=t';
t=t-1977;
near5year= [2011:1:2015]';
near5year=near5year-1977;
%阻滯增長增長模型
beta0 = [566 0.22 12567];
yzuzhi =@(beta,t) beta(3)./(1+(beta(3)./beta(1)-1).*exp(-beta(2)*t));
beta=nlinfit(t,x,yzuzhi,beta0)
xzuzhi=yzuzhi(beta,near5year);
%馬爾薩斯增長模型
ymuthus=@(r,t) r(1)*exp(r(2)*t);
r0=[5.7,0.1];
a=nlinfit(t,x,ymuthus,r0)
xmuthus=ymuthus(a,near5year);
subplot(1,2,1);
plot([t;near5year],[x;xmuthus]);
title('馬爾薩斯增長模型');
subplot(1,2,2);
plot([t;near5year],[x;xzuzhi]);
title('阻滯增長增長模型');
%預(yù)測未來五年
disp(['馬爾薩斯增長模型預(yù)測未來五年工資是',num2str(xmuthus')]);
disp(['阻滯增長增長模型預(yù)測未來五年工資是',num2str(xzuzhi')]);
clear all;clc
x=[566 632 745 755 769 789 985 1110 ...
1313 1428 1782 1920 2150 2292 2601 3149 ...
4338 5145 5809 6271 6854 7656 8772 10007 ...
11374 12567 14332 16614 19228 22844 26404 29638 ...
32074]';
n = 33;
t = (1:n)';
beta0 = [5.3 0.22 400]; %[x0, r, xm]
[beta, R, J] = nlinfit(t,x,'logisfun',beta0); %非線性擬合
py = beta(3)./(1+(beta(3)./beta(1)-1).*exp(-beta(2)*t)); %預(yù)測各年工資
p5 = beta(3)./(1+(beta(3)./beta(1)-1).*exp(-beta(2)*(n:n+5)')); %預(yù)測五年后工資
subplot(2,1,1)
plot(1:n, x, '*', 1:n, py,n:n+5,p5,'o'); %對比圖
title('阻滯增長增長模型');
xx = x(1:n);
t=[ones(n,1),(1:n)'];
y=log(xx(1:n));
[b,bint,r,rint,stats] = regress(y,t)
x0 = exp(b(1)) %參數(shù)x0
r = b(2)%參數(shù)r
py_ = x0*exp(r*(1:n+5)'); %預(yù)測數(shù)據(jù)
subplot(2,1,2)
plot(1:n, xx, 'k+', 1:n+5, py_); %對比圖
title('馬爾薩斯增長模型');
8.練習(xí)（傳染病模型）
假設(shè)某一種傳染病正在流行，現(xiàn)在希望建立適當(dāng)?shù)臄?shù)學(xué)模型，利用已經(jīng)掌握的一些數(shù)據(jù)資料對該傳染病進(jìn)行有效地研究，用于對其傳播蔓延進(jìn)行必要的控制，減少人民生命財(cái)產(chǎn)損失?，F(xiàn)假設(shè)某地區(qū)人口分布約為5萬人，且有約為5千人已經(jīng)感染并積極治療，每個病人每天接觸約為5人，且每天治愈率約為150%。考慮如下問題，建立適當(dāng)?shù)臄?shù)學(xué)模型：
- 該傳染病不具有免疫性，試估算病情穩(wěn)定時對應(yīng)的天數(shù)，并繪制出病人與健康者人數(shù)隨時間變化的分布曲線；
- 該傳染病具有免疫性，試估算病情穩(wěn)定時對應(yīng)的天數(shù)，并繪制出病人，健康者與免疫者人數(shù)隨時間變化的分布曲線.
8.1 傳染病無免疫性
病人治愈成為健康人，健康人可以再次被感染（SIS模型）
假設(shè)：
- 總?cè)藬?shù)N不變，病人和健康人比例分別為：i(t),s(t)
- 每個病人每天有效接觸人數(shù)為 λ \lambda λ,且使接觸的健康人致病
- 病人每天治愈的比例為 μ \mu μ
建模：
N
[
i
(
+
Δ
t
)
?
i
(
t
)
]
=
λ
s
(
t
)
N
i
(
t
)
Δ
t
?
μ
N
i
(
t
)
Δ
t
N[i(+\Delta t)-i(t)] = \lambda s(t)Ni(t)\Delta t -\mu Ni(t)\Delta t
N[i(+Δt)?i(t)]=λs(t)Ni(t)Δt?μNi(t)Δt
得到：
{
d
i
d
t
=
λ
i
s
?
u
i
i
(
0
)
=
i
0
\left\{\begin{aligned} \cfrac{di}{dt}&=\lambda is-ui \\ \\ i(0)&=i_0\\ \end{aligned}\right.
??????????dtdi?i(0)?=λis?ui=i0??
clear all;
clc;
xspan = [0 5];
y0 = 0.1;
[x,y] = ode45(@(x,y) 5*y*(1-y)-1.5*y, xspan, y0);
y1 = 1-y(:)
plot(x,y,'r',x,y1,'b')
由圖可以看出，需要兩天病情才能穩(wěn)定。
8.2 傳染病有免疫性
病人治愈后即移除感染系統(tǒng)，稱為移出者（SIR模型）
假設(shè)：
- 總?cè)藬?shù)N不變, 病人、健康人、移出者的比例分別為 i ( t ) , s ( t ) , r ( t ) i(t),s(t),r(t) i(t),s(t),r(t)
- 病人每天有效接觸人數(shù)為 λ \lambda λ，治愈的比例為 μ \mu μ,接觸數(shù)為 σ = λ / μ \sigma = \lambda / \mu σ=λ/μ
建模：
s
(
t
)
+
i
(
t
)
+
r
(
t
)
=
1
s(t)+i(t)+r(t)=1
s(t)+i(t)+r(t)=1
需要建立
s
(
t
)
,
i
(
t
)
,
r
(
t
)
的
兩
個
方
程
s(t),i(t),r(t)的兩個方程
s(t),i(t),r(t)的兩個方程：
N
[
i
(
t
+
Δ
t
)
?
i
(
t
)
]
=
λ
s
(
t
)
N
i
(
t
)
Δ
t
?
μ
N
i
(
t
)
Δ
t
N
[
s
(
t
+
Δ
t
)
?
s
(
t
)
]
=
?
λ
s
(
t
)
N
i
(
t
)
Δ
t
N[i(t+\Delta t)-i(t)] = \lambda s(t)Ni(t)\Delta t -\mu Ni(t)\Delta t \\ \\ N[s(t+ \Delta t)-s(t)] =-\lambda s(t)Ni(t)\Delta t\\
N[i(t+Δt)?i(t)]=λs(t)Ni(t)Δt?μNi(t)ΔtN[s(t+Δt)?s(t)]=?λs(t)Ni(t)Δt
{
d
i
d
t
=
λ
i
s
?
u
i
d
s
d
t
=
?
λ
s
i
i
(
0
)
=
i
0
,
s
(
0
)
=
s
0
\left\{\begin{aligned} \cfrac{di}{dt}&=\lambda is-ui\\ \cfrac{ds}{dt}&=-\lambda si \\ i(0)&=i_0,s(0)=s_0\\ \end{aligned}\right.
??????????????dtdi?dtds?i(0)?=λis?ui=?λsi=i0?,s(0)=s0??
無法求出s(t),i(t)的解析解，利用matlab求出數(shù)值解文章來源：http://www.zghlxwxcb.cn/news/detail-779897.html
function y = infect(t,x)
lamp = 5;
u = 1.5;
y = [lamp*x(1)*x(2)-u*x(1); -lamp*x(1)*x(2)];
end
x0 = [0.9, 0.1]';
[t,x] = ode45('infect', [0,4], x0);
%調(diào)用變步長4階5級Runge-Kutta-Felhberg法計(jì)算
y = 1-x(:,1)-x(:,2)
plot(t,x(:,1),'r', t,x(:,2),'b',t,y,'k');
grid on
由圖可以看出，大概4天后，傳染病趨于穩(wěn)定。文章來源地址http://www.zghlxwxcb.cn/news/detail-779897.html
到了這里，關(guān)于【數(shù)學(xué)建模\MATLAB】掌握用Matlab求解微分方程問題的文章就介紹完了。如果您還想了解更多內(nèi)容，請?jiān)谟疑辖撬阉鱐OY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關(guān)文章，希望大家以后多多支持TOY模板網(wǎng)！
