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

3年前浏览2556

过冷水诚挚邀请你加入Matlab仿真秀官方交流群进行Matlab学习、问题咨询、 Matlab相关资料下载，群号：927550334

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

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

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

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.

我们再来看一个案例：

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

%第一种解微分方程的方法

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);