方程式求根

方程式解析解

想要获得方程式的解析解, 首先要能够明确表示方程的解析式. 表示方程式可以通过符号变量来实现, 符号变量是一种新接触的变量, 可以使用syms或者sym来定义.

1
2
3
syms x % 定义符号变量 x
%% 等价于
x=sym('x') % 定义符号变量 x

成功定义后, 可以在工作区中发现一个新的变量x, 其类型为 sym

image-20200714210440014

中学数学教我们x+x+x应该等于3x, 此处, 如果定义x是符号变量的话, 在matlab中输入x+x+x 也会得到一个3*x的解析解.

1
2
3
syms x
x + x + x
(x + x + x)/4

image-20200714210807755

其运算结果ans的类型也是符号变量. 也就是说, 由符号变量导出的解析式的类型也是符号变量.因此对于解析式y=x2−2x−8y=x^2-2x-8可以通过先定义符号变量x, 再输入方程关系式得到y.

1
2
syms x
y = x^2-2*x-8

下面正式介绍如何求方程式的解析解.

中学数学告诉我们, 求方程式的根与求函数零点问题是等价的, 因此一切方程式求根问题, 都可以转化为f(x)=0f(x)=0的问题. matlab中使用solve()函数来求解函数的零点或者叫方程式的根. solve()需要两个参数, 第一个参数输入函数解析式, 第二个参数指定哪个符号变量是自变量, 下面给出一些使用示例.

求y=x⋅sin⁡(x)−x=0y=x \cdot \sin (x)-x=0 的根.

1
2
3
4
5
6
7
%% 第一种方式
syms x
solve('x*sin(x)-x', x)
%% 第二种方式
syms x
y = x*sin(x)-x;
solve(y, x)

solve()函数也可以用于方程组求解.

{x−2y=5x+y=6\left\{\begin{array}{c} x-2 y=5 \\ x+y=6 \end{array}\right.

1
2
3
4
syms x y            % 定义x 和 y 两个变量
eq1 = x - 2*y - 5;  % 分别写出两个解析式
eq2 = x + y - 6;
A = solve(eq1,eq2,x,y) % 前两个参数先给出解析式, 再给出哪些是自变量

上述代码执行后, 在命令行窗口中给出如下答案:

image-20200714212541614

我们可以在工作区中查看A的类型为结构体类型, 包括x 和y两个字段.

image-20200714212533459

可以通过执行A.x和A.y来查看解的具体数值.

image-20200714212745952

solve函数不仅可以用于求解具体数值, 也可以用于求公式解.

求ax2−b=0ax^2-b=0的公式解.

1
2
syms x a b
solve('a*x^2-b')

image-20200714212357259

需要指出的是, 上面的代码中, solve函数并没有指定哪个是未知数, 哪些是已知数. 在不指定哪个变量是未知数的情况下, 系统将默认x 是未知数. 如果我们把问题改成用a 和x 来表示b, 此时就需要在solve函数中指出哪个是未知数.代码应该改为:

1
2
syms x a b
solve('a*x^2-b','b')

image-20200714213040731

练习:

  1. 使用其他变量表示x : (x−a)2+(y−b)2=r2(x-a)^{2}+(y-b)^{2}=r^{2}
  2. 求矩阵[abcd]\left[\begin{array}{ll} a & b \\ c & d \end{array}\right]的逆矩阵.

参考代码:

1
2
3
%% 第一题
syms x a y b r
solve('(x-a)^2+(y-b)^2-r^2','x')

运行结果:

image-20200714213638304

1
2
3
4
5
6
7
syms a b c d A B C D
eq1=a*A+b*C-1;
eq2=a*B+b*D;
eq3=c*A+d*C;
eq4=c*B+d*D-1;
Ans=solve(eq1,eq2,eq3,eq4,'A','B','C','D');
[Ans.A Ans.B;Ans.C Ans.D]

运行结果:

image-20200714214236947

上一篇文章中学习了求微积分的数值解, 这里介绍一下求微积分解析式的函数. diff在上一篇文章中介绍了其在处理向量时的作用, 在处理符号变量时, diff将直接给出函数的解析式的导函数解析式.

求y=4x5y=4x^5的导函数:

1
2
3
syms x
y = 4*x^5;
yprime = diff(y)

image-20200714214819096

相对应的不定积分的求法使用int()函数.

求z=∫ydx=∫x2exdx,z(0)=0z=\int y d x=\int x^{2} e^{x} d x, \quad z(0)=0

1
2
3
4
syms x; y = x^2*exp(x);
z = int(y);
z = z-subs(z, x, 0) % subs函数实现的对解析式z中的x替换成0 也就是说subs(z,x,0)会得到函数z(x)的常数项的值.
% 因此, z减去其常数项的值 就满足了z(0)=0的条件. 将其赋值给z即可.

image-20200714215449608

练习:

  1. 计算导数.

    f(x)=ex2x3−x+3,dfdx=?f(x)=\frac{e^{x^{2}}}{x^{3}-x+3}, \quad \frac{d f}{d x}=?

    g(x)=x2+xy−1y3+x+3,∂f∂x=?g(x)=\frac{x^{2}+x y-1}{y^{3}+x+3}, \quad \frac{\partial f}{\partial x}=?

  2. 计算定积分.

    ∫010x2−x+1x+3dx\int_{0}^{10} \frac{x^{2}-x+1}{x+3} d x

参考代码:

1
2
3
4
5
6
7
8
%% 计算f'(x)
syms x
f=(exp(x^2))/(x^3-x+3);
diff(f)
%% 计算g'(x)
syms x y
g=(x^2+x*y-1)/(y^3+x+3);
diff(g,'x') % 指出求导对象
1
2
3
4
%% 计算定积分 此处利用牛顿-莱布尼兹公式
syms x
y=int((x^2-x+1)/(x+3));
subs(y,x,10)-subs(y,x,0)

积分结果为:

image-20200714220922978

方程式数值解

在数值分析那门课程中, 介绍了求解方程的数值解的一些方法, 数值分析是一门很难的课程, 以后有机会可能会尝试学习一下.

在中学的学习中, 我们所接触的方程 包括一元方程组, 一元二次方程等等都可以导出其求根公式, 也就是说那些方程的根与其未知数的系数有关系, 通过后来的学习, 也认识到一元三次方程也存在公式解. 那么, 对于一元四次方程, 五次方程等等高次方程有没有公式解呢? 经过复杂的证明, 5次方程以上的方程没有公式解! 对于第一次听到这个消息的小伙伴可能也会像我一样有些理解不了, 方程式的根竟然不能用系数来表示? 但是事实就是这样. 因此, 研究方程的数值解法就显得很有必要.

方程式求根的数值解法常见的有:

  • 二分法
  • 牛顿法

二分法思想比较简单, 不做过多介绍. 这里直接给出算法流程图.

image-20200714222744551

下面介绍一下牛顿法的逼近过程是怎样的.

牛顿法首先需要从某一点x0x_0出发, 计算f(x0)f(x_0)以及在点f(x0)f(x_0)处的切线.

记此切线与x轴相交于点x1x_1, 计算f(x1)f(x_1)以及在点f(x1)f(x_1)处的切线,

记此切线与x轴相交于点x2x_2, . . . .

随着迭代次数的增加, f(xn)f(x_n)的数值会逐渐减小 , 当∣f(xn)∣|f(x_n)|小于精度时, 即可停止迭代.其算法示意图如下:

image-20200714222643308

算法流程图:

image-20200714222703434

下表中给出二分法与牛顿法的主要区别:

image-20200714222859261


matlab中求方程的数值解通常使用fzero()或者fsolve(), 这两个函数的传入参数都是函数句柄

fzero()函数相当于使用二分法求解, 因此函数必须穿过x轴, 对于y=x2y=x^2这样的函数, fzero()无法求解, 因为区间端点的函数值的乘积大于0.

fsolve()相当于使用牛顿法求解, 因此除了要给出求解函数之外, 还需要给出迭代的初始点. 下面分别举例做介绍.

求f(x)=1.2x+0.3+x⋅sin⁡(x)f(x)=1.2 x+0.3+x \cdot \sin (x) 的零点

1
2
f2 = @(x) (1.2*x+0.3+x*sin(x));
fsolve(f2,0) % f2为函数句柄 0是初始迭代点

运行结果:

image-20200714223858675

这两个函数还有一些可选参数.

1
2
3
4
f=@(x)x.^2
options=optimset('MaxIter',1e3,'TolFun',1e-10); % 指定迭代步数 和 精度
fsolve(f,0.1,options)
fzero(f,0.1,options)

运行结果:

image-20200714224148739

对于多项式, 可以使用roots函数来求方程式的数值解, 例如:

求解f(x)=x5−3.5x4+2.75x3+2.125x2−3.875x+1.25f(x)=x^{5}-3.5 x^{4}+2.75 x^{3}+2.125 x^{2}-3.875 x+1.25

1
roots([1 -3.5 2.75 2.125 -3.875 1.25])

运行结果:

image-20200714224423997

练习:

求方程组f(x,y)={2x−y−e−x−x+2y−e−yf(x, y)=\left\{\begin{array}{c} 2 x-y-e^{-x} \\ -x+2 y-e^{-y} \end{array}\right.的根, 使用初始迭代点为(x,y)=(−5,−5)(x, y)=(-5,-5)

参考代码:

1
2
3
4
F = @(x) [2*x(1) - x(2) - exp(-x(1));
         -x(1) + 2*x(2) - exp(-x(2))];
x0 = [-5;-5];
[x,fval] = fsolve(F,x0)

运行结果:

image-20200714230628731

下面比较一下数值计算与解析计算的各自的优势与缺点.

image-20200714230807876

上面介绍的数值计算的方法都是基于迭代的, 还有一类问题处理方式与上面的方法恰恰相反, 那就是递归问题. 递归函数就是指在函数体内又调用了函数本身.

比如求解n!n!时 , 就可以采用递归算法进行求解.

1
2
3
4
5
6
7
8
9
function output = fact(n)
% fact recursively finds n!
if n==1
	output = 1;
else
	output = n * fact(n-1);
end
end

资料链接

参考视频 参考讲义