利用 Maxima 求解常微分方程

本文最初写于 2010-10-16 于 sohu 博客，这次博客搬家一起搬到这里来。

版权所有，转载请注明出处。

 

含带导数符号或带微分符号的未知函数的方程称为微分方程。
如果在微分方程中未知函数是一个变元的函数，这样的微分方程称为常微分方程。

1 一阶、二阶常微分方程的通解

Maxima 可以求解很多种类的常微分方程。
对于可以给出闭式解的一阶和二阶常微分方程，Maxima 会试图求出其精确解。
下面给出三个简单的例子。

(%i1) eq1:'diff(y,x)+y=x;
sol1:ode2(eq1,y,x);

 

 

(%i3) eq2:'diff(y,x,2)+y=sin(3*x);
sol2:ode2(eq2,y,x);

(%i5) eq3:'diff(y,x,3)+y=0;
sol3:ode2(eq3,y,x);

上面的例子用了ode2函数来求解常微分方程。
在定义方程时，微分函数diff之前有一个单引号（‘），这表示让Maxima只给出形式上的输出，并不真的进行计算。
这是因为我们这里只要列出方程，并不想让Maxima真的求导。
sol1 中的%c 和 sol2 中的 %k1 %k2 是任意常数。
ode2函数只能求解一阶和二阶常微分方程，第三个例子给出的是一个三阶常微分方程，无法求解，因此输出 false。

 

2 初值问题

函数ic1 (solution, xval, yval)和ic2 (solution, xval, yval, dval)分别用来解一阶和二阶微分方程的初值问题，其中solution是用ode2解得的通解，xval和yval分别是自变量和因变
量的初值，dval是因变量一阶导数的初值。

(%i7) ic1(sol1,x=0,y=1);

(%i8) ic2(sol2,x=0,y=1,'diff(y,x)=-1);

3 边值问题

函数bc2 (solution, xval_1, yval_1, xval_2, yval_2)用来求解二阶微分方程的边值问题，
其中solution是ode2解得的通解，xval_1、yval_1、xval_2和yval_2分别为自变量和因变量在第一点和第二点的取值。

(%i9) eq4:'diff(y,x,2)+diff(y,x)/2+3*y=0;
sol4:ode2(eq4,y,x);
bc2(sol4,x=0,y=1,x=5,y=-1);

 

4 利用Laplace变换法求解常微分方程(组)

如果待求解的常微分方程(组)是线性常系数的。则可以利用Laplace变换法来求解。
Maxima 中也提供了相应的求解函数 desolve()，desolve()函数既可以求解ODE 方程，也可以求解ODE方程组。函数的基本形式如下。
desolve (eqn, y)
desolve ([eqn_1, ..., eqn_n], [y_1, ..., y_n])
这里待解函数不能只写变量名（例如y），而需要明确写出对自变量的依赖关系（例如y(x)）。

下面是一个简单的例子：

(%i12) kill(all);
eq5:'diff(y(x),x,3)+y(x)=0;
desolve(eq5,y(x));

如果初值是已知的，可以使用atvalue()命令来提供初值。
如果提供了足够的初值条件，再用的desolve()函数求解时积分常数自然就可以确定了。

(%i3) atvalue(y(x),x=0,1);
desolve(eq5,y(x));
atvalue('diff(y(x),x),x=0,-1);
desolve(eq5,y(x));
atvalue('diff(y(x),x,2),x=0,0);
desolve(eq5,y(x));

下面给出一个常微分方程组求解的例子。

(%i9) kill(all);
eq6:'diff(f(x),x,2)+'diff(g(x),x)+3*f(x)=15*exp(-x);
eq7:'diff(g(x),x,2)-4*'diff(f(x),x)+3*g(x)=16*sin(x);
ode_syst:[eq6,eq7];
atvalue(f(x),x=0,35);
atvalue('diff(f(x),x),x=0,-48);
atvalue(g(x),x=0,27);
atvalue('diff(g(x),x),x=0,-55);
desolve(ode_syst,[f(x),g(x)]);

下面是试验部分。

(%i9) kill(all);
eq1:'diff(y(x),x)+y(x)=x;
sol1:desolve(eq1,y(x));

(%i3) atvalue(y(x),x=1,1);
desolve(eq1,y(x));

(%i5) ic1(y=(%c+1)*exp(-x)+x-1,x=1,y(x)=1);

说明 desolve 函数提供的初值必须是x=0 处的。
ic1 不能用来直接求解 desolve 函数的结果。必须要人为的处理一下结果的形式。这一点上确实不方便。