B站台湾大学郭彦甫｜MATLAB 学习笔记｜11 方程式求根 Root_Finding

MATLAB学习笔记（11 方程式求根 Root_Finding）

1. 提出问题

假设现在有一个数学函数 $f(x)=x_2-2x-8$，需要找到一个点 $x_0$，使得 $f(x_0)=0$，使用 MATLAB 来解决这个问题：

• 求函数的解析解
• 使用图像来表达
• 数值方法求解

2. 找符号形式的根

2.1 使用 solve() 找符号形式的根 (Symbolic approach)

syms x	% 创建symbols的意思
y = x*sin(x)-x;
solve(y, x) % slove(y,x)中y表示equation,x表示variable，默认等式y=0

>> Class11
 
ans =
 
   0
pi/2

2.2 求方程组

用符号方式求下列方程组：

$$\begin{cases} x-2y=5\\ x+y=6\\ \end{cases}$$

syms x y
eq1 = x - 2*y - 5;
eq2 = x + y - 6;
A = solve(eq1,eq2,x,y)

exercise (P7)

题目 1：

使用符号方式求解下列等式对 $x$ 的解：

$$\left( x-a \right) ^2+\left( y-b \right) ^2=r^2$$

解 1：

%% 题目1
clear
clc
syms x y a b r
f = (x-a)^2 + (y-b)^2 - r^2;
solve(f,x) 

ans =
 
a + (b + r - y)^(1/2)*(r - b + y)^(1/2)
a - (b + r - y)^(1/2)*(r - b + y)^(1/2)

题目 2：

使用符号方式求解矩阵的逆：

$$A=\left[ \begin{matrix} a& b\\ c& d\\ \end{matrix} \right]$$

解 2：

%% 题目2
clear
clc
syms a b c d 
A = [a b; c d];
inv(A)  % 对A求逆

ans =
 
[ d/(a*d - b*c), -b/(a*d - b*c)]
[-c/(a*d - b*c),  a/(a*d - b*c)]

2.3 求差分：diff() （得到符号形式结果）

example:

计算下列函数的导数：

$$y=4x^5$$

syms x
y = 4*x^5;
yprime = diff(y)	%yprime 表示y的导数

exercise (P8)

题目：

$$\begin{align} f\left( x \right) &=\frac{e^{x^2}}{x^3-x+3}, \ \ \ \frac{df}{dx} =? \\ g\left( x \right) &=\frac{x^2+xy+1}{y^3+x+3}, \ \ \frac{\partial g}{\partial x}=?\text{（PPT上输错了）} \end{align}$$

解：

clear
clc
syms x y
f = (exp(x^2))/(x^3-x+3);
g = (x^2+x*y+1)/(y^3+x+3);
diff(f,x)
diff(g,x) 

ans =
 
(2*x*exp(x^2))/(x^3 - x + 3) - (exp(x^2)*(3*x^2 - 1))/(x^3 - x + 3)^2
 
 
ans =
 
(2*x + y)/(y^3 + x + 3) - (x^2 + y*x + 1)/(y^3 + x + 3)^2

2.4 符号积分

example:

计算以下函数的积分：

$$z=\int{ydx=\int{x^2e^xdx, \ \ \text{且满足} z\left( 0 \right) =0}}$$

对 $y$ 直接积分得到 $z=e^x(x^2-2x+2)$，此时的 $z(0)=2$，为使得 $z(0)=0$，因此（我们所需要的） $z$ = （现在求得的） $z$ - 2

syms x; 
y = x^2*exp(x);
z = int(y); 
z = z-subs(z, x, 0);
% subs(s,old,new)返回s的副本，其中所有出现的old将会被new取代。

exercise (P9)

题目：

计算：

$$\int_0^{10}{\frac{x^2-x+1}{x+3}dx}$$

解：

clear
clc
syms x 
y = (x^2-x+1)/(x+3);
F = int(y,x,0,10)

 
F =
 
log(302875106592253/1594323) + 10

3. 找数值形式的根 (Numerical approach)

3.1 使用 fsolve() 求根

example:

计算：

$$f\left( x \right) =1.2x+0.3+x\cdot \sin \left( x \right)$$

f2 = @(x) (1.2*x+0.3+x*sin(x));
fsolve(f2,0)    % f2是句柄，0是计算初始值

exercise (P14)

题目：

找到下面方程组的根：

$$f\left( x,y \right) =\begin{cases} 2x-y-e^{-x}\\ -x+2y-e^{-y}\\ \end{cases}$$

初始值 $(x,y)=(-5,-5)$

解：

clear
clc

fun = @eqas;
x0 = [-5 -5];
answer = fsolve(fun,x0)

function F = eqas(x0) 
% 在eqas的输入为1个数组时才能计算出来，还需要注意F
x = x0(1); y = x0(2);
F(1) = 2*x-y-exp(-x);
F(2) = -x+2*y-exp(-y);
end


answer =

    0.5671    0.5671

3.2 使用 fzero()

example：

计算：$y=x^2$ 的根

%方法1
f=@(x)x.^2；
fzero(f,0.1)	%fzero()必须保证函数经过x轴，才能算出来

fsolve(f,0);	%作为对比

ans =

   NaN
   
ans =

     0

%方法2
f=@(x)x.^2
options=optimset('MaxIter',1e3,'TolFun',1e-10);
%1e3表示迭代次数（1000），1e-10表示误差
fsolve(f,0.1,options)
fzero(f,0.1,options)

f =

  包含以下值的 function_handle:

    @(x)x.^2


Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

<stopping criteria details>

ans =

   1.9532e-04

正在退出 fzero: 将终止搜索包含符号变化的区间
 因为在搜索期间遇到 NaN 或 Inf 函数值。
(-1.37296e+154 处的函数值为 Inf。)
请检查函数或使用其他起始值重试。

ans =

   NaN

3.3 求多项式的根: roots()

example：

找到下面多项式的根：

$$f\left( x \right) =x^5-3.5x^4+2.75x^3+2.125x^2-3.875x+1.25$$

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

ans =

   2.0000 + 0.0000i
  -1.0000 + 0.0000i
   1.0000 + 0.5000i
   1.0000 - 0.5000i
   0.5000 + 0.0000i

3.4 数值寻根方法进阶

以数值的方式寻根, 其方法可分为两种主要类型:

1. Bracketing methods (交叉法)

   从一个包含根的区间开始, 例如 bisection method (二分法)

2. Open methods (开型法)

   从一个或多个初始预测点开始, 例如 Newton-Raphson method (牛顿-拉普森方法)

寻根将会迭代进行，直到满足某些标准：精度或者迭代次数

3.4.1 二分法 (Bisection method)

• 假设条件

  • $f(x)$ 在 $[l,u]$ 上连续
  • $f\left( l \right) \cdot f\left( u \right) <0$

• 算法

将空间分成一半，判断两段端点的函数值的积的正负，选择函数积为负的一半空间继续使用二分法

Loop

1. $r=\left( l+u \right) /2$（找中点）

2. 如果 $f\left( r \right) \cdot f\left( u \right) <0$，确定新的间隔 $[r,u]$

   如果 $f\left( l \right) \cdot f\left( u \right) <0$，确定新的间隔 $[l,r]$

End

• 二分法算法流程图

3.4.2 牛顿-拉普森方法 (Newton-Raphson method)

• 假设条件

  • $f(x)$ 连续
  • $f'(x)$ 已知

• 算法

根据已知的一个点的 $x_0,f(x_0),f'(x_0)$ 找到斜率的延长线与 X 轴的交点 $x_1$，作为新的计算点，以此进行迭代计算

$$\begin{aligned} \frac{f\left( x_0 \right)}{x_0-x_1}&=f'\left( x_0 \right) \\ \Longrightarrow x_1&=x_0-\frac{f\left( x_0 \right)}{f'\left( x_0 \right)}\ \ \ \text{(根据$x_0,f(x_0),f'(x_0)$ 来求$x_1$)} \end{aligned}$$

Loop

$\ \ \ x_{n+1}=x_n-\frac{f\left( x_n \right)}{f'\left( x_n \right)}$

End

• 算法流程图

4. 递归方法求根

4.1 介绍递归函数

函数可以自己调用自己

example:

函数如下：（以阶乘递归函数为例）

$$\begin{aligned} n ! &=n \times(n-1) ! \\ &=n \times(n-1) \times(n-2) ! \\ &=n \times(n-1) \times(n-2) \times(n-3) ! \\ &=n \times(n-1) \times(n-2) \times \cdots \end{aligned}$$

function output = fact(n)
% fact recursively finds n!
if n==1
output = 1;
else
output = n * fact(n-1);
end
end

该函数包括一个递归情况和一个基数情况，当函数到达基数情况时停止