MATLAB求高阶微分方程数值解

注:MATLAB不能直接求高阶微分方程的数值解，需要变量替换，将其转化为一阶微分方程组。

例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

再调用

%例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%解析解求导，即为y'
[x,y1]=ode45('wf1',[-3,4],[2,4]);%数值解
%注意: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的数值解和解析解')
legend('y数值解','y解析解')
subplot(2,1,2)
plot(x,y1(:,2),'r*',x,y3,'b+')
title('y的一阶导数的数值解和解析解')
legend('y的一阶导数的数值解','y的一阶导数的解析解')

例2

该一阶微分方程组的二阶微分方程形式为:

y''=1000(1-y²)y'-y

首先，编写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

再调用

%例2
%刚性方程组
%使用 ode45 以及默认的相对和绝对误差容限（分别为 1e-3 和 1e-6）解算此方程组特别慢，需要好几分钟才能完成解算和绘制。
[t,y] = ode15s('wf2',[0 3000],[2 0]);
plot(t,y(:,1),'-o')

MATLAB求高阶微分方程数值解.rar