带你用matlab轻松搞定微分方程

之前过冷水有和大家分享热传导方程求解的方法，其本质上是微分方程的问题。考虑大多数读者对微分方程求解方法比较陌生，所以过冷水本期简单普及一下微分方程的求解问题。

关于微分方程你需要了解：含有未知的函数及其某些阶的导数以及其自变量本身的方程称为微分方程。如果未知函数是一元函数，则称为常微分方程。如果未知函数是多元函数，则称为偏微分方程。联系一些未知函数的一组微分方程称为微分方程组。微分方程中出现的未知函数的导数的最高阶称为微分方程的阶。

有些微分方程比较简单可直接通过积分求解。例如一阶常系数线性常微分方程：

syms y(x) a b
eqn = diff(y,x) == a*y+b;
S = dsolve(eqn)
S =
-(b - C3*exp(a*x))/a

syms y(x) a b
eqn = diff(y,x,2)+2*diff(y,x,1)+4*y ==0；
S = dsolve(eqn)
eqn(x) =
4*y(x) + 2*diff(y(x), x) + diff(y(x), x, x) == 0
S =
C4*exp(-x)*cos(3^(1/2)*x) + C5*exp(-x)*sin(3^(1/2)*x)

演示了两个比较简单的微分方程用符号解微分方程的方法解出通解，在我们实际问题中少数特殊方程可用初等积分法求解外，大部分微分方程无显示解，应用主要依靠数值解法。考虑一阶常微分方程组初值问题：

其中y=(y₁,y₂,...,y_m)^T，f=(f₁,f₂,...,f_m)^T，y₀=(y₁₀,y₂₀,...,y_m0)^T，所谓数值解，就是寻求解y(t)在一些列离散节点t0<t1<t2<...<t_n<t_f  上的近似值y_k(k=0,1,...,n).称h_k=t_k+1-t_k为步长，已知：

求其数值解。自己根据差分方程思想编程如下：

clear all
warning off
feature jit off
f=inline('y-2*x/y','x','y');
a=0;b=1;h=0.1;
n=(b-a)./h;
x=zeros(1,n+1);
y=zeros(1,n+1);
y(1)=1;
for i=1:n+1
    x(i)=a+(i-1)*h;
    y(i+1)=y(i)+h*f(x(i),y(i));
end
y=y(1:n+1);

一般来讲符号法的运算会比单纯的数值运算可具有科学准确性。因为该问题比较简单，可以采用符号微分法求解，用符号计算为对比看差分法数值运算精度如何。代码如下：

%符号法解微分方程
y1=dsolve('Dy=y-2*x/y','y(0)=1','x')
%绘图对比
figure1 = figure;
axes1 = axes('Parent',figure1);
hold(axes1,'on');
plot(x,y,'DisplayName','差分法','LineWidth',2);
h1=ezplot(y1,[0,1]);
set(h1,'DisplayName','符号法','LineWidth',2);
xlabel('$x$','FontWeight','bold','Interpreter','latex');
ylabel('$f(x)$','Interpreter','latex');
title('差分法&符号法 解微分方程对比');
xlim(axes1,[0 1]);
ylim(axes1,[0.93232679283255 1.79972401473633]);
set(axes1,'FontSize',14,'FontWeight','bold');
legend1 = legend(axes1,'show');
set(legend1,'Position',[0.148120235966763 0.82820510022613 0.127868850472195 0.0933333308498066]);

我们再来看一个案例：

拟采用两种符号运算方法，两种数值运算方法。代码如下：

%第一种解微分方程的方法
syms y(x) a b
eqn1 = diff(y,x,1)==-2*y(x)+2*x^2+2*x;
eqn2 = y(0)==1;
y1= dsolve([eqn1,eqn2]);
%第二种解微分方程的方法
y2=dsolve('Dy=-2*y+2*x^2+2*x','y(0)=1','x');
%第三种解微分方程的方法
warning off
feature jit off
f=inline('-2*y+2*x^2+2*x','x','y');
a=0;b=0.5;h=0.0001;
n=(b-a)./h;
x=zeros(1,n+1);
y=zeros(1,n+1);
y(1)=1;
for i=1:n+1;
    x3(i)=a+(i-1)*h;
    y(i+1)=y(i)+h*f(x(i),y(i));
end
y3=y(1:n+1);
%第四种微分方程数值解
[x4,y4]=ode23(f,[0 0.5],1);
%绘图
figure1 = figure;
axes1 = axes('Parent',figure1);
hold(axes1,'on');
h1=ezplot(y1,[0,0.5]);h2=ezplot(y2,[0,0.5]);
set(h1,'DisplayName','$y_1 =exp(-2x) + x^2$','LineWidth',2);
set(h2,'DisplayName','$y_2=2x + 2exp(-2x) - x^2 - 1$','LineWidth',2);
plot(x3,y3,'DisplayName','$(x_3,y_3)$','LineWidth',2);
plot(x4,y4,'DisplayName','$y_4=ode23(f(x))$','MarkerFaceColor',[0.494117647409439 0.184313729405403 0.556862771511078],'MarkerSize',8,'Marker','o','LineStyle','none');
xlabel('$x$','FontSize',16,'Interpreter','latex');
ylabel('$y_4$','FontSize',20,'Interpreter','latex');
title('微分方程不同解法数值解对比图');
xlim(axes1,[0 0.5]);
ylim(axes1,[0.43978341778928 1.0459754645536]);
set(axes1,'FontSize',14,'LineWidth',2);
legend1 = legend(axes1,'show');
set(legend1,'Position',[0.604420982252212 0.73205553519439 0.220018931648527 0.166277798138944],'Interpreter','latex','FontSize',14);

本期推文过冷水就是想讲一下简单的微分方程求解的方法，让大家足够解决常见问题就OK了！至于那些复杂问题，万丈高楼平地起，Monte Carlo算法不也讲了好几期的吗？敬请期待下期的复杂偏微分方程组的求解方法。

往期回顾>>>>>>

欢迎各路英雄豪杰来搞

积分变量替换到legendre微分变换

数值计算——MATLAB数值积分原理详讲

数值优化—三种复杂函数数值积分方法实例演示

非线性方程组求解迭代算法&图像寻初始值讲解

互动专区

在matlab爱好者公众号中，回复“QQ”加入公众号专属Q群；回复“原创”获取小编原创代码；回复“星球”加入资源分享园地知识星球。

如需转载，请在公众号中回复“转载”获取授权，未经授权擅自搬运抄袭的，必将追究其责任！