实例讲解Matlab求解二阶常系数微分方程

今天，又有群友问了一个关于二阶常系数微分方程得求解问题，一共两个二阶常系数微分方程，还是个分段函数，要用Matlab求解该方程。相当于，一个二阶常系数微分方程，不同自变量范围得到两个方程。每个分端函数，都有两个边界条件，一共四个边界条件。

这里我们要用到得Matlab函数分别有syms变量定义、diff微分方程、dsolve微分方程求解等。下面，正式开始求解得过程。

问题描述如下：

古哥的方法

我们分析，这种方程是由两个相同自变量x不同范围组成得两个分段函数，需要分开边界条件，求解得到两个函数。最后，再用两个边界条件联立，求解两个函数中未确认得常数变量（Matlab中，求解结果如果有不确定得结果，会尝试用C来表示）。我们在Matlab m文件中，输入如下代码：

clc; 
clear all;
syms u1(x) u2(x) beta1 beta2 r l q e
eqn1=diff(u1,x,2)-diff(u1,x).*beta1==r.*q.*(l./2-x);
eqn2=diff(u2,x,2)-diff(u2,x).*beta2==r.*q.*(l./2-x);
du1=diff(u1,x);
du2=diff(u2,x);
cond1=[u1(0)==0];
cond2=[u2(l/2)==0];
u1=dsolve(eqn1,cond1)
u2=dsolve(eqn2,cond2)

运行结果如下：

 
u1 =
 
C3 - C3*exp(beta1*x) + (q*r*x^2)/(2*beta1) - (q*r*x*(beta1*l - 2))/(2*beta1^2)
 
 
u2 =
 
C6 - (exp(-(beta2*l)/2)*exp(beta2*x)*(8*C6*beta2^2 - q*r*beta2*l^2 + 4*q*r*l))/(8*beta2^2) + (q*r*x^2)/(2*beta2) - (q*r*x*(beta2*l - 2))/(2*beta2^2)
 
>> 

我们可以看到，表达式中，有已定义常量：beta1、q、r、l、beta2，以及已定义变量x，以及未定义常数C3、C6。即，我们再用最后两个边界条件，代入u1、u2得表达式，就可以确认C3、C6的值。

由于，对方没有给我beta1、q、r、l、beta2这几个常量的具体值，这里我就不方便继续求解确认C3、C6的值了。根据两个条件确认两个变量（C3、C6），可以确认该方程可以得到确认的解。因此，我们判断问题得到解决。

别人的方法

后续来了，我一般都晚上有时间，对方下午的时候又付费让别人写了一段程序来求解，看起来写法更加简单，代码如下：

clear
clc
syms beta1 beta2 gama q L e 
y=dsolve('D2u1-beta1^2*u1=gama*q*(L/2-x)','D2u2-beta2^2*u2=gama*q*(L/2-x)+beta2^2*u2','Du1(0)=0','Du2(L/2)=0','Du1(e)=Du2(e)','u1(e)=u2(e)','x');
u1=simplify(y.u1)
u2=simplify(y.u2)

运行结果如下：

 
u1 =
 
(exp(beta1*x)*((exp(beta1*e)*(4*beta2^3*gama*q*exp(-e*(beta1 - 2*2^(1/2)*beta2)) - 4*beta2^3*gama*q*exp(2^(1/2)*L*beta2 - beta1*e) + 2^(1/2)*beta1^3*gama*q*exp(2^(1/2)*L*beta2) + 2^(1/2)*beta1^3*gama*q*exp(2*2^(1/2)*beta2*e) - 2*2^(1/2)*beta1^3*gama*q*exp((2^(1/2)*beta2*(L + 2*e))/2) + 2*2^(1/2)*beta1*beta2^2*gama*q*exp(2^(1/2)*L*beta2 - beta1*e) + 2*L*beta1*beta2^3*gama*q*exp(2^(1/2)*L*beta2) - L*beta1^3*beta2*gama*q*exp(2^(1/2)*L*beta2) - 2*2^(1/2)*beta1*beta2^2*gama*q*exp(2^(1/2)*L*beta2) - 4*beta1*beta2^3*e*gama*q*exp(2^(1/2)*L*beta2) + 2*beta1^3*beta2*e*gama*q*exp(2^(1/2)*L*beta2) - 2*L*beta1*beta2^3*gama*q*exp(2*2^(1/2)*beta2*e) + L*beta1^3*beta2*gama*q*exp(2*2^(1/2)*beta2*e) - 2*2^(1/2)*beta1*beta2^2*gama*q*exp(2*2^(1/2)*beta2*e) + 4*beta1*beta2^3*e*gama*q*exp(2*2^(1/2)*beta2*e) - 2*beta1^3*beta2*e*gama*q*exp(2*2^(1/2)*beta2*e) + 2*2^(1/2)*beta1*beta2^2*gama*q*exp(-e*(beta1 - 2*2^(1/2)*beta2))))/(2*beta1^2*beta2^2*(2*beta2*exp(2*beta1*e + 2^(1/2)*L*beta2) + 2*beta2*exp(2^(1/2)*L*beta2) - 2*beta2*exp(2*2^(1/2)*beta2*e) - 2*beta2*exp(2*e*(beta1 + 2^(1/2)*beta2)) - 2^(1/2)*beta1*exp(2*2^(1/2)*beta2*e) + 2^(1/2)*beta1*exp(2*e*(beta1 + 2^(1/2)*beta2)) + 2^(1/2)*beta1*exp(2*beta1*e + 2^(1/2)*L*beta2) - 2^(1/2)*beta1*exp(2^(1/2)*L*beta2))) + (gama*q*exp(-beta1*x)*(2*beta1*x - L*beta1 + 2))/(4*beta1^2)))/beta1 - (exp(-beta1*x)*((gama*q*exp(beta1*x)*(L*beta1 - 2*beta1*x + 2))/(4*beta1^2) - (gama*q*exp(beta1*e)*(4*beta2^3*exp(beta1*e + 2^(1/2)*L*beta2) - 4*beta2^3*exp(e*(beta1 + 2*2^(1/2)*beta2)) + 2^(1/2)*beta1^3*exp(2^(1/2)*L*beta2) + 2^(1/2)*beta1^3*exp(2*2^(1/2)*beta2*e) - 2*2^(1/2)*beta1^3*exp((2^(1/2)*beta2*(L + 2*e))/2) + 2*2^(1/2)*beta1*beta2^2*exp(beta1*e + 2^(1/2)*L*beta2) + 2*L*beta1*beta2^3*exp(2^(1/2)*L*beta2) - L*beta1^3*beta2*exp(2^(1/2)*L*beta2) - 2*2^(1/2)*beta1*beta2^2*exp(2^(1/2)*L*beta2) - 4*beta1*beta2^3*e*exp(2^(1/2)*L*beta2) + 2*beta1^3*beta2*e*exp(2^(1/2)*L*beta2) - 2*L*beta1*beta2^3*exp(2*2^(1/2)*beta2*e) + L*beta1^3*beta2*exp(2*2^(1/2)*beta2*e) - 2*2^(1/2)*beta1*beta2^2*exp(2*2^(1/2)*beta2*e) + 4*beta1*beta2^3*e*exp(2*2^(1/2)*beta2*e) - 2*beta1^3*beta2*e*exp(2*2^(1/2)*beta2*e) + 2*2^(1/2)*beta1*beta2^2*exp(e*(beta1 + 2*2^(1/2)*beta2))))/(2*beta1^2*beta2^2*(2*beta2*exp(2*beta1*e + 2^(1/2)*L*beta2) + 2*beta2*exp(2^(1/2)*L*beta2) - 2*beta2*exp(2*2^(1/2)*beta2*e) - 2*beta2*exp(2*e*(beta1 + 2^(1/2)*beta2)) - 2^(1/2)*beta1*exp(2*2^(1/2)*beta2*e) + 2^(1/2)*beta1*exp(2*e*(beta1 + 2^(1/2)*beta2)) + 2^(1/2)*beta1*exp(2*beta1*e + 2^(1/2)*L*beta2) - 2^(1/2)*beta1*exp(2^(1/2)*L*beta2)))))/beta1
 
 
u2 =
 
(2^(1/2)*exp(2^(1/2)*beta2*x)*((gama*q*exp(-2^(1/2)*beta2*x)*(2*2^(1/2)*beta2*x - 2^(1/2)*L*beta2 + 2))/(8*beta2^2) - (gama*q*exp(2^(1/2)*beta2*e)*(4*beta2^3*exp(2*beta1*e) - 8*beta2^3*exp(beta1*e) - 2*beta1^2*beta2 + 4*beta2^3 - 2*beta1^2*beta2*exp(2*beta1*e) + 2^(1/2)*beta1^3*exp(2*beta1*e + (2^(1/2)*L*beta2)/2 - 2^(1/2)*beta2*e) + 2*beta1^2*beta2*exp(2*beta1*e + (2^(1/2)*L*beta2)/2 - 2^(1/2)*beta2*e) - 2^(1/2)*beta1^3*exp((2^(1/2)*beta2*(L - 2*e))/2) + 2*beta1^2*beta2*exp((2^(1/2)*beta2*(L - 2*e))/2) - 2*L*beta1*beta2^3 + L*beta1^3*beta2 + 4*beta1*beta2^3*e - 2*beta1^3*beta2*e + 2*L*beta1*beta2^3*exp(2*beta1*e) - L*beta1^3*beta2*exp(2*beta1*e) - 4*beta1*beta2^3*e*exp(2*beta1*e) + 2*beta1^3*beta2*e*exp(2*beta1*e)))/(2*beta1^2*beta2^2*(2*beta2*exp(2*beta1*e + 2^(1/2)*L*beta2) + 2*beta2*exp(2^(1/2)*L*beta2) - 2*beta2*exp(2*2^(1/2)*beta2*e) - 2*beta2*exp(2*e*(beta1 + 2^(1/2)*beta2)) - 2^(1/2)*beta1*exp(2*2^(1/2)*beta2*e) + 2^(1/2)*beta1*exp(2*e*(beta1 + 2^(1/2)*beta2)) + 2^(1/2)*beta1*exp(2*beta1*e + 2^(1/2)*L*beta2) - 2^(1/2)*beta1*exp(2^(1/2)*L*beta2)))))/(2*beta2) - (2^(1/2)*exp(-2^(1/2)*beta2*x)*((gama*q*exp(2^(1/2)*beta2*x)*(2^(1/2)*L*beta2 - 2*2^(1/2)*beta2*x + 2))/(8*beta2^2) + (gama*q*exp(2^(1/2)*beta2*e)*(4*beta2^3*exp(2*beta1*e + 2^(1/2)*L*beta2) - 8*beta2^3*exp(beta1*e + 2^(1/2)*L*beta2) + 4*beta2^3*exp(2^(1/2)*L*beta2) - 2*beta1^2*beta2*exp(2^(1/2)*L*beta2) - 2^(1/2)*beta1^3*exp(2*beta1*e + (2^(1/2)*L*beta2)/2 + 2^(1/2)*beta2*e) + 2*beta1^2*beta2*exp(2*beta1*e + (2^(1/2)*L*beta2)/2 + 2^(1/2)*beta2*e) + 2^(1/2)*beta1^3*exp((2^(1/2)*beta2*(L + 2*e))/2) + 2*beta1^2*beta2*exp((2^(1/2)*beta2*(L + 2*e))/2) - 2*beta1^2*beta2*exp(2*beta1*e + 2^(1/2)*L*beta2) + 2*L*beta1*beta2^3*exp(2*beta1*e + 2^(1/2)*L*beta2) - L*beta1^3*beta2*exp(2*beta1*e + 2^(1/2)*L*beta2) - 4*beta1*beta2^3*e*exp(2*beta1*e + 2^(1/2)*L*beta2) + 2*beta1^3*beta2*e*exp(2*beta1*e + 2^(1/2)*L*beta2) - 2*L*beta1*beta2^3*exp(2^(1/2)*L*beta2) + L*beta1^3*beta2*exp(2^(1/2)*L*beta2) + 4*beta1*beta2^3*e*exp(2^(1/2)*L*beta2) - 2*beta1^3*beta2*e*exp(2^(1/2)*L*beta2)))/(2*beta1^2*beta2^2*(2*beta2*exp(2*beta1*e + 2^(1/2)*L*beta2) + 2*beta2*exp(2^(1/2)*L*beta2) - 2*beta2*exp(2*2^(1/2)*beta2*e) - 2*beta2*exp(2*e*(beta1 + 2^(1/2)*beta2)) - 2^(1/2)*beta1*exp(2*2^(1/2)*beta2*e) + 2^(1/2)*beta1*exp(2*e*(beta1 + 2^(1/2)*beta2)) + 2^(1/2)*beta1*exp(2*beta1*e + 2^(1/2)*L*beta2) - 2^(1/2)*beta1*exp(2^(1/2)*L*beta2)))))/(2*beta2)
 
>> 

他这种方式，直接将所有边界条件，一起用了，结果有点长。