MATLAB求高階微分方程數(shù)值解
e小白網(wǎng)址:www.e-xiaobai.com
注:MATLAB不能直接求高階微分方程的數(shù)值解，需要變量替換，將其轉(zhuǎn)化為一階微分方程組。
例1
首先，要將二階微分方程降為一階微分方程組，如下圖所示:
首先，編寫wf1.m文件
function dy=wf1(x,y)
dy=zeros(2,1);
dy(1)=y(2);
dy(2)=2*x/(1+x^2)*y(2);
end再調(diào)用
%例1
%dsolve('(1+x^2)*D2y=2*x*Dy','y(-3)=2,Dy(-3)=4','x')%求得解析解
%(2*x*(x^2 + 3))/15 + 34/5%解析解，即為y
%(2*x^2)/5 + 2/5%解析解求導(dǎo)，即為y'
[x,y1]=ode45('wf1',[-3,4],[2,4]);%數(shù)值解
%注意:y1變量中有兩列值，第1列為y1的值，即為題目中y的值。
%第2列為y2的值，即為題目中y'的值。
y2=(2.*x.*(x.^2 + 3))/15 + 34/5;
y3=(2.*x.^2)/5 + 2/5;
figure(1)
subplot(2,1,1)
plot(x,y1(:,1),'r*',x,y2,'b+')
title('y的數(shù)值解和解析解')
legend('y數(shù)值解','y解析解')
subplot(2,1,2)
plot(x,y1(:,2),'r*',x,y3,'b+')
title('y的一階導(dǎo)數(shù)的數(shù)值解和解析解')
legend('y的一階導(dǎo)數(shù)的數(shù)值解','y的一階導(dǎo)數(shù)的解析解')
例2
該一階微分方程組的二階微分方程形式為:
首先，編寫wf2.m文件
function dydx=wf2(t,y)
dydx=zeros(2,1);
dydx(1)=y(2);
dydx(2)=1000*(1-y(1)^2)*y(2)-y(1);
end再調(diào)用
%例2
%剛性方程組
%使用 ode45 以及默認(rèn)的相對(duì)和絕對(duì)誤差容限（分別為 1e-3 和 1e-6）解算此方程組特別慢，需要好幾分鐘才能完成解算和繪制。
[t,y] = ode15s('wf2',[0 3000],[2 0]);
plot(t,y(:,1),'-o')
標(biāo)簽：
