4.8  微分方程
微分方程是数值计算中常见的问题，MATLAB提供了多种函数来计算微分方程的解。
4.8.1  常微分方程
众所周知，对一些典型的常微分方程，能求解出它们的一般表达式，并用初始条件确定表达式中的任意常数。但实际中存在有这种解析解的常微分方程的范围十分狭窄，往往只局限在线性常系数微分方程(含方程组)，以及少数的线性变系数方程。对于更加广泛的、非线性的一般的常微分方程，通常不存在初等函数解析解。由于实际问题求解的需要，求近似的数值解成为了解决问题的主要手段。常见的求数值解的方法有欧拉折线法、阿当姆斯法、龙格-库塔法与吉尔法等。其中由于龙格-库塔法的精度较高，计算量适中，所以使用的较广泛。
数值解的最大优点是不受方程类型的限制，即可以求任何形式常微分方程的特解(在解存在的情况下)，但是求出的解只能是数值解。
1．龙格-库塔方法简介
对于一阶常微分方程的初值问题，在求解未知函数时，在点的值是已知的，并且根据高等数学中的中值定理，应有：
一般而言，在任意点，有：
当确定后，根据上述递推公式能算出未知函数在点的一列数值解：
当然在递推的过程中同样存在着一个误差累计的问题，实际计算中的递推公式一般都进行过改造，龙格-库塔公式为：
其中：
2．龙格-库塔法的实现
基于龙格-库塔法，MATLAB提供了ode系列函数求常微分方程的数值解。常用的有ode23 和ode45函数，其调用语法如下。
(1)[t,y]=ode23(filename,tspan,y0)：采用了二阶、三阶龙格-库塔法进行计算。
(2)[t,y]=ode45(filename,tspan,y0)：采用了四阶、五阶龙格-库塔法进行计算。
其中filename是定义f(t,y)的函数文件名，该函数文件必须返回一个列向量。tspan的形式为[t0,tf]，表示求解区间。y0是初始状态列向量。t和y分别给出时间向量和相应的状态向量。
这两个函数分别采用了二阶、三阶龙格-库塔法和四阶、五阶龙格-库塔法，并采用自适应变步长的求解方法，即当解的变化较慢时采用较大的步长，从而使得计算的速度很快；当解的变化较快时步长会自动地变小，从而使得计算的精度很高。
【例4-45】  设有初值问题：
试求其数值解，并和精确解相比较，精确解为()。
首先要建立微分方程所对应的函数文件myodefun.m，文件内容如下：
function y=myodefun(t,y)                      
%  建立函数文件myodefun.m
y=(y^2-t-2)/(4*(t+1));
建立myodefun函数之后，就可以调用ode23函数求解微分方程。
>> t0=0;
>> tf=10;
>> y0=2;
>> [t,y]=ode23 ('myodefun',[t0,tf],y0);       %  求数值解
>> y1=sqrt(t+1)+1;                              %  求精确解
>> plot(t,y,'k.',t,y1,'r')
通过图形比较，数值解用黑色圆点表示，精确解用红色实线表示，如图4-13所示。
【例4-46】  求下面无劲度系统微分方程组的数值解。
为了求解方程，首先要建立方程的m文件。本例中不妨建立名为rigid.m的函数文件，此文件用以描述给出的方程组，文件的内容如下：
function dy = rigid(t,y)
dy = zeros(3,1);                 %  一个列向量
dy(1) = y(2) * y(3);
dy(2) = -y(1) * y(3);
dy(3) = -0.51 * y(1) * y(2);
本例中，我们通过odeset函数对误差进行控制，另外在时间[0 12]进行求解，0时刻初始条件向量为[0 1 1]。
>> options =odeset('RelTol',1e-4,'AbsTol',[1e-4 1e-4 1e-5]);%  误差控制
>> [T,Y] = ode45(@rigid,[0 12],[0 11],options);               %  求数值解
>>plot(T,Y(:,1),'-',T,Y(:,2),'-.',T,Y(:,3),'.')               %  绘制结果图
得到的结果如图4-14所示。
  
图4-13  常微分方程结果图 
图4-14  常微分方程数值

【作者声明】
　　本文所有文字均为作者原创，所有图片均为作者本人亲自拍摄或制作。
　　版权所有，仅供阅读欣赏，禁止任何单位或个人以任何形式对本文的文字或图片进行包括但不限于复制、转载、引用、抄袭、截图、模仿、翻译、印刷等之中的一项或多项的行为。禁止将本文用于商业用途。
　　作者保留所有权利，请尊重作者的劳动成果，谢谢合作。
前言
　　我们知道，fx-991CN X计算器的“方程/函数”模式、SOLVE功能都可以用来求解多项式方程，前者可以求解标准形式的一元二~四次方程，后者使用牛顿迭代法（切线法）求解。但是这两个功能都有局限性，就是不能求解系数为复数的方程。而且在fx-991CN X的复数模式中，SOLVE功能是无法使用的，因此我们必须想办法利用现有的功能来求解。
数学基础与功能预备知识
牛顿迭代公式： 
　　在计算模式下，SOLVE功能就是使用这个公式求解的，我们只需要输入方程，再指定一个初始值，计算器就能够按照初始值来计算方程的一个根。例如求解方程 
 ，我们只需要按[菜单]、[1]进入计算模式，输入方程，按[SHIFT]、[CALC]进入SOLVE功能，指定一个迭代初始值（例如5）按[=]确认，再次按下[=]即可得到一个解“x=3”。
如果不是很清楚SOLVE功能的使用，具体可以参考以前的这篇文章：
電卓院亜紀良：fx-991CN X方程（组）求解功能的使用​zhuanlan.zhihu.com
复数模式：按[菜单]、[2]进入复数模式，可以进行复数的四则运算、整数次幂、共轭复数、实部/虚部提取、代数形式(a+bi)-极坐标形式(r∠θ)相互转换等等。
如果不太了解复数功能，具体可以参考以前的这篇文章：
電卓院亜紀良：使用科学计算器计算复数与相量（基础篇）​zhuanlan.zhihu.com
在复数模式下求解复系数多项式方程
　　在复数模式中，虽然SOLVE功能不能使用，但CALC功能仍然可以使用，那么我们可以使用CALC功能输入迭代式，间接求解复系数的多项式方程。
【例】求解复系数方程 
【解】根据牛顿迭代公式，我们需要计算出函数 
的导数 
 ，多项式函数的导数是很容易计算的。这里给的是一个标准形式的多项式方程，如果遇到的多项式方程不是标准形式，最好先经过手算变形整理成标准形式。那么 
 的导数就是
所以迭代公式就是
在计算器上就是输入
按[菜单]、[2]进入复数模式，输入上述迭代式，注意等号按[ALPHA]、[CALC]输入红色的等号，然后按[CALC]。
这里需要输入迭代的初始值，一般输入一个实数即可，注意不要让原方程对应的函数的导数值为0，否则会出现数学错误。这里就以0为初始值，按[=]即可开始迭代；如果想要使用其他值作为初始值，那么就再输入一个数，按[=]确认之后，再按[=]开始迭代。迭代的时候会交替显示迭代当前值与赋值界面。
按多次[=]后，即可发现x的值逐渐趋于不变，这时候即可认为已经迭代到了一个根：
按[S-D]可以将实部和虚部显示为两行：
如果要得到其他的根，在迭代完成的时候继续按[=]进入赋值界面，重新指定一个初始值，例如这里指定的是 
 ：
这样，这个复系数一元二次方程就解出来了：
后记
　　这篇文章书写起因是一个提问，有关复系数方程求解的问题不止一次出现过，虽然带有CAS的图形计算器能够直接求解，但在只能使用fx-991CN X这样的计算器的考试中，利用现有的工具，加以简单的手算，搭配一些算法，就可以轻松地解决计算器不能直接求解的问题。另外还要感谢 @朽木家的优姬 、 @小林露露 的讨论与支持。

使用MATLAB求解方程求根——学习笔记

碎碎念：终于参加完了某比赛，连续大约摸了两天的鱼，就在昨天由于自己的操作失误，亲手将电脑给烧了，这就是上天在暗示我是时候加油为接下来的两场比赛和一个考试努力啦~下面浅谈一下MATLAB各种求解方程的根的方法。

求方程的根

法一：

求方程的全部根。

    在MATLAB命令窗口输入：

    p=[1,7,0,9,-20];   %建立多项式系数向量

    x=roots(p)         %求根

法二：

使用solve函数：

s=solve(f,v)：求方程关于指定自变量的解；

s=solve(f)：求方程关于默认自变量的解。

其中 f可以是用字符串表示的方程，或符号表达式；若f 中不含等号，则表示解 方程 f=0。

解方程 x^3-3*x+1=0

>>syms x; f=x^3-3*x+1;

>>s=solve(f,x)

>>s=solve(‘x^3-3*x+1’,‘x’)

>>s=solve(‘x^3-3*x+1=0’,‘x’)

Solve函数还可以用来求解方程组！

solve(f1,f2,...,fN,v1,v2,...,vN)

求解由 f1,f2,...,fN 确定的方程组关于 v1, v2,...,vN 的解。

解方程组：

>>[x,y,z]=solve(‘x+2*y-z=27’,‘x+z=3’, ...

                ‘x^2+3*y^2=28’,‘x’,‘y’,‘z’)

注意：输出变量的顺序要书写正确！

Solve在得不到解析解时，会给出数值解。

非线性方程组求解：

fzero(f,x0)：求方程f=0在x0附近的根。

几点说明：

方程可能有多个根，但是fzero只能给出距离x0最近的一个根；
	
若x0是一个标量，则fzero先找出一个包含x0的区间，使得f在这个区间两个端点上的值异号，然后再在这个区间内寻找方程f=0的根；如果找不到这样的区间，则返回NaN。
	
若x0是一个2维向量，则表示在[x0(1),x0(2)]区间内求方程的根，此时必须满足f在这两个端点上的值异号。
	
由于fzero是根据函数是否船业横轴来决定零点。因此它无法确定函数曲线仅触及横轴但是不穿越的零点，如|sin(x)|的所有零点。
	
函数中的f是一个函数句柄，通过以下方式给出：
字符串形式：fzero(‘x^3-3*x+1’,2);

通过@调用的函数句柄：fzero(@sin,4);

f不能用符号表达式！
 

线性方程组求解：

linsolve(A,b)：解线性方程组    

 

>> A=[1 2 –1; 1 0 1; 1 3 0];

>> b=[2;3;8];

>> X=linsolve(A,b)

注意这里的b是列向量。

求解方程函数小结

roots(p)：多项式的所有零点，p是多项式系数向量。

fzero(f,x0)：求f=0在x0附近的根，f是函数句柄，可以由字符串给出或使用@， 但不能是符号表达式！

solve(f,v)：求方程关于指定自变量的解，f可以是用字符串表示的方程，或符号表 达式，若不含等号表示f=0；

也可解方程组(包含非线性)；得不到解析解时，给出数值解。

linsolve(A,b)：解线性方程组。

在matlab中解方程组是很方便的
例如，对于代数方程组Ax=b(A为系数矩阵，非奇异)的求解，MATLAB中有两种方法：

(1)x=inv(A)*b — 采用求逆运算解方程组；

(2)x=A\b — 采用左除运算解方程组。

例:

x1+2x2=8 

2x1+3x2=13

>>A=[1,2;2,3];b=[8;13];

>>x=inv(A)*b 

x = 

   2.00 

   3.00 

>>x=A\b

x = 

  2.00

  3.00；

即二元一次方程组的解x1和x2分别是2和3。


对于同学问到的用matlab解多次的方程组，有符号解法，方法是：先解出符号解,然后用vpa(F,n)求出n位有效数字的数值解.具体步骤如下：

第一步:定义变量syms x y z ...;

第二步:求解[x,y,z,...]=solve('eqn1','eqn2',...,'eqnN','var1','var2',...'varN');

第三步:求出n位有效数字的数值解x=vpa(x,n);y=vpa(y,n);z=vpa(z,n);...。

如：解二（多）元二（高）次方程组：

x^2+3*y+1=0

y^2+4*x+1=0

解法如下：

>>syms x y;

>>[x,y]=solve('x^2+3*y+1=0','y^2+4*x+1=0');

>>x=vpa(x,4);

>>y=vpa(y,4);

结果是：

x = 

    1.635+3.029*i

    1.635-3.029*i

    -.283

   -2.987

y = 

    1.834-3.301*i

    1.834+3.301*i

    -.3600

   -3.307。

二元二次方程组，共4个实数根；



还有的同学问，如何用matlab解高次方程组（非符号方程组）？举个例子好吗？

解答如下：

基本方法是：solve(s1,s2,…,sn,v1,v2,…,vn)，即求表达式s1,s2,…,sn组成的方程组，求解变量分别v1,v2,…,vn。

具体例子如下：

x^2 + x*y + y = 3

x^2 - 4*x + 3 = 0

解法：

>> [x,y] = solve('x^2 + x*y + y = 3','x^2 - 4*x + 3 = 0')

运行结果为 

x =

     1 3

y =

     1 -3/2


即x等于1和3；y等于1和-1.5


或

>>[x,y] = solve('x^2 + x*y + y = 3','x^2 - 4*x + 3= 0','x','y')

 x =

     1  3

 y =

     1 -3/2

结果一样，二元二方程都是4个实根。


通过这三个例子可以看出，用matlab解各类方程组都是可以的，方法也有多种，只是用到解方程组的函数，注意正确书写参数就可以了，非常方便。
用Matlab解方程组的时候，发现它不能自动代入系数的值。比如

说如下的程序;

a=4;

x=solve('a*x=4',x);

怎么解决？下面是解决办法，很简单。
法1：
symsa x b;

b = solve('a*x=4', x);

a = 4;

b=eval(b);
法2：
symsa x b;a = 4;

b = solve('a*x=4', x);
b=subs(b)