《数值计算方法》系列总目录
第一章 误差序列实验
第二章 非线性方程f(x)=0求根的数值方法
第三章 CAD模型旋转和AX=B的数值方法
第四章 插值与多项式逼近的数值计算方法
第五章 曲线拟合的数值方法
第六章 数值微分计算方法
第七章 数值积分计算方法
第八章 数值优化方法
第二章
一、算法原理
1、波尔查诺二分法
求解非线性方程的根的时候，二分法主要适用于在已知方程的根的大致区间的情况下，通过二分法可将区间内的端点逐步逼近零点，直到得到一个任意小的包含零点的间隔。二分法判定过程的第1步是选择中点。分析可能存在的三种情况。
1. List item
2. 如果f(b)和f©符号相反，则在区间[c,b]内存在零点。
3. 如果f©=0，则c是零点。
如果情况1或情况2，则表示找到了一个比
lim
n
→
∞
c
n
=
r
\mathop {\lim }\limits_{n \to \infty } {c_n} = r
n→∞limcn=r原先区间范围小一半的区间，它包含根并称之为对区间进行压缩。为了持续此过程，需要对新的更小区间进行重新标号。重复执行直到区间足够小，在误差允许的范围内，我们便认为找到了方程根的近似值。
设
f
∈
C
[
a
,
b
]
f \in C[a,b]
f∈C[a,b]，且存在数
r
∈
[
a
,
b
]
r \in [a,b]
r∈[a,b]满足
f
(
r
)
=
0
f(r) = 0
f(r)=0，如果f(a)和f(b)的符号相反，且
{
c
n
}
n
=
0
∞
\{ {c_n}\} ^\infty _{n = 0}
{cn}n=0∞表示二分法生成的中点序列，则
∣
r
−
c
n
∣
≤
b
−
a
2
n
+
1
|r - {c_n}| \le \frac{{b - a}}{{{2^{n + 1}}}}
∣r−cn∣≤2n+1b−a,其中n=0,1…。这样序列
{
c
n
}
n
=
0
∞
\{ {c_n}\} ^\infty _{n = 0}
{cn}n=0∞收敛到零点x=r即可表示为
lim
n
→
∞
c
n
=
r
\mathop {\lim }\limits_{n \to \infty } {c_n} = r
n→∞limcn=r。
通过迭代求解区间中点的值，重新确定区间，当函数值的绝对值小于给定误差范围的时候，找到了近似解。
2、试值法
试值法与二分法的原理基本相同，但是收敛速度相对较快。与二分法条件一样，找到给定区间的割线 与X轴的交点，把这一交点作为近似值。迭代过程中，比较函数值的正负，并重新确定区间。通过迭代，最终在给定误差范围内得到方程的近似解。
通过对斜率不同的计算方式，可以得到每次迭代中的近似值的表达形式，与二分法的判别条件相同。
- 如果 f ( a ) f(a) f(a)和 f ( c ) f(c) f(c)符号相反，则在区间[a,c]内存在零点。
- 如果 f ( b ) f(b) f(b)和 f ( c ) f(c) f(c)符号相反，则在区间[c,b]内存在零点。
- 如果 f ( c ) = 0 f(c)=0 f(c)=0，则c是零点。
最终，通过比较 ∣ f ( c ) ∣ |f(c)| ∣f(c)∣的与给定误差范围的大小，确定迭代的终止条件，并认为找到了可以接受的方程的近似解。
3、牛顿拉夫森法
牛顿拉夫森法比二分法和试值法收敛速度更快。基本原理是假设初始值
p
0
{p_0}
p0在根
p
p
p附近。函数
y
=
f
(
x
)
y=f(x)
y=f(x)的图形与x轴相交于点
(
p
,
0
)
(p,0)
(p,0)。而且点
(
p
0
,
f
(
p
0
)
)
({p_0},f({p_0}))
(p0,f(p0))位于靠近点
(
p
,
0
)
(p,0)
(p,0)的曲线上。将
p
1
{p_1}
p1定义为曲线在点
(
p
0
,
f
(
p
0
)
)
({p_0},f({p_0}))
(p0,f(p0))的切线与x轴的交点,通过图像的显示则可以看到
p
1
{p_1}
p1比
p
0
{p_0}
p0更靠近
p
p
p。如果写出切线L的两种表达式。可得到
p
1
{p_1}
p1和
p
0
{p_0}
p0相关的方程。
m
=
0
−
f
(
p
0
)
p
1
−
p
0
m
=
f
′
(
p
0
)
{\rm{ }}m = \frac{{0 - f({p_0})}}{{{p_1} - {p_0}}}{\rm{ }}m = f'({p_0})
m=p1−p00−f(p0)m=f′(p0)
p
1
=
p
0
−
f
(
p
0
)
f
′
(
p
0
)
{\rm{ }}{p_1} = {p_0} - \frac{{f({p_0})}}{{f'({p_0})}}
p1=p0−f′(p0)f(p0)
则
p
k
=
g
(
p
k
−
1
)
=
p
k
−
1
−
f
(
p
k
−
1
)
f
′
(
p
k
−
1
)
{\rm{ }}{p_k} = g({p_{k - 1}}) = {p_{k - 1}} - \frac{{f({p_{k - 1}})}}{{f'({p_{k - 1}})}}
pk=g(pk−1)=pk−1−f′(pk−1)f(pk−1)
二、实验内容及核心算法代码
1、波尔查诺二分法
利用二分法求解函数的零点，通过二分法的实验原理，我们设置起始区间值为 =0, =2。方程的根位于这个区间内。通过不断的取区间的中点 ，并求出中点的函数值f( )，并与区间两端的函数值进行比较，进而来重新确定区间的中点。通过比较函数值的绝对值| f( )|与设置误差的大小，从而确定方程最终收敛的近似值。
其中，程序的核心算法代码为
float GetfunrootbyBinsearch(float (*f)(float), float x0, float x1, int& nIter, const int MaxIter, const float eps, float** valset,bool& ifOuepreci)
{
float val0, val1, xh, valh;
nIter = 0;//the raw index
//if the initial is the root,return it
val0 = f(x0);
if (fabs(val0) < eps)
{
return x0;
}
val1 = f(x1);
if (fabs(val1) < eps)
{
return x1;
}
//if set the error initial,program exit
if (f(x0) * f(x1) > 0)
{
cout << "set the error initial!program exit!" << endl;
exit(0);
}
//Porchano binary search method
do
{
xh = (x0 + x1) / 2;
valh = f(xh);
valset[nIter][0] = nIter + 1; valset[nIter][1] = x0; valset[nIter][2] = xh; valset[nIter][3] = x1; valset[nIter][4] = valh;
if (val0 * valh > 0)
{
x0 = xh;
val0 = valh;
}
else
{
x1 = xh;
val1 = valh;
}
++nIter;//move the index
if (nIter == MaxIter)
{
cout << "cannot find the root within the MaxIter num!" << endl;
ifOuepreci = true;
break;
}
} while (fabs(valh) > eps);
return xh;
#endif
这段代码的作用是实现，不断求取与判断的区间中点的函数值与给定误差的大小，从而作为判断循环的条件。并将每一步得到数据保存在一个二维指针指向的空间内。当所求得的函数值小于给定误差时，作为取得方程的近似解，并返回这一近似解。
2、试值法求解xsin(x)-1=0原理实现
与二分法求解方程的近似解的条件一致，设置起始区间值为
a
0
a_0
a0=0,
b
0
b_0
b0=2，利用试值法的原理的表达式
c
n
=
b
n
−
f
(
b
n
)
(
b
n
−
a
n
)
f
(
b
n
)
−
f
(
a
n
)
{c_n} = {b_n} - \frac{{f({b_n})({b_n} - {a_n})}}{{f({b_n}) - f({a_n})}}
cn=bn−f(bn)−f(an)f(bn)(bn−an)
进行迭代，从而不断生成新的取值区间，从而缩小近似解的取值区间。判断的所取区间的近似解的函数值与给定误差的大小，从而求得可以接受的方程的近似解。
其核心代码如下
float GetCk(float (*f)(float), float x0, float x1)
{
float meval1, meval2;
meval1 = f(x1) * (x1 - x0);
meval2 = f(x1) - f(x0);
return x1 - meval1 / meval2;
}
float GetfunrootbyFSP(float (*f)(float), float x0, float x1, int& nIter, const int MaxIter, const float eps, float** valset, bool& ifOuepreci)
{
float val0, val1, xh, valh;
nIter = 0;//the raw index
//if the initial is the root,return it
val0 = f(x0);
if (fabs(val0) < eps)
{
return x0;
}
val1 = f(x1);
if (fabs(val1) < eps)
{
return x1;
}
//if set the error initial,program exit
if (f(x0) * f(x1) > 0)
{
cout << "set the error initial!program exit!" << endl;
exit(0);
}
//Porchano binary search method
do
{
xh = GetCk(f,x0,x1);
valh = f(xh);
valset[nIter][0] = nIter + 1; valset[nIter][1] = x0; valset[nIter][2] = xh; valset[nIter][3] = x1; valset[nIter][4] = valh;
if (val0 * valh > 0)
{
x0 = xh;
val0 = valh;
}
else
{
x1 = xh;
val1 = valh;
}
++nIter;//move the index
if (nIter == MaxIter)
{
cout << "cannot find the root within the MaxIter num!" << endl;
ifOuepreci = true;
break;
}
} while (fabs(valh) > eps);
return xh;
}
3、牛顿拉弗森法求解sqrt(5)的近似值原理实现
利用牛顿方法原理中推导出来的表达式，通过求得所给函数在A等于5时候的零点，从而求出根号5的近似值。迭代函数的简要推导如下所示。与之前面的二分法和试值法迭代过程相似。通过简化迭代函数的表达式，从而减少迭代过程中计算机运算的次数，从而减小误差，同时通过传递引用的方式，在一个封装的函数中同时得到多个相关的结果。
f
(
x
)
=
x
2
−
A
f(x) = {x^2} - A
f(x)=x2−A
g
(
x
)
=
x
−
f
(
x
)
f
′
(
x
)
=
x
−
x
2
−
A
2
x
g(x) = x - \frac{{f(x)}}{{f'(x)}} = x - \frac{{{x^2} - A}}{{2x}}
g(x)=x−f′(x)f(x)=x−2xx2−A
g
(
x
)
=
x
+
A
/
x
2
g(x) = \frac{{x + A/x}}{2}
g(x)=2x+A/x
p
k
=
p
k
−
1
+
A
p
k
−
1
2
{\rm{ }}{p_k} = \frac{{{p_{k - 1}} + \frac{A}{{{p_{k - 1}}}}}}{2}
pk=2pk−1+pk−1A
其核心代码如下
float GetSqrtrootbyNewtonMethod(float (*g)(float, float), const float eps, float* pn, float p0, int& nIter, float A, int nMaxIter, bool& ifFindRoot)
{
float err;
pn[0] = p0;
nIter = 0;
if (fabs(pn[0] - g(A, pn[0])) < eps)
{
return pn[0];
}
do
{
++nIter;
pn[nIter] = g(A, pn[nIter - 1]);
err = (pn[nIter] - pn[nIter - 1]) / pn[nIter];
if (nIter >= nMaxIter)
{
ifFindRoot = false;
cout << "cannot find the root within the precision!" << endl;
break;
}
} while (fabs(err) > eps);
return pn[nIter];
4、牛顿拉弗森法求解炮弹飞行轨迹问题原理实现
通过分析炮弹在飞行过程中的运动状态，并考虑在飞行过程中的所受的阻力，得到炮弹在飞行过程中的运动学方程，且为一个非线性的方程。通过对炮弹进行在水平方向和竖直方向上的运动分析，利用Newton-Raphson方法，从而计算出炮弹从发出到落地所需要的时间，即f(t)=0的近似解，从而得出炮弹飞行过程中的飞行时间，最高高度，水平位移距离等多个参数。通过分析可得如下表达式。
y
=
υ
y
t
−
16
t
2
,
a
n
d
x
=
υ
x
{\rm{ }}y = {\upsilon _y}t - 16{t^2},{\rm{ }}and{\rm{ }}x = {\upsilon _x}
y=υyt−16t2,andx=υx
y
=
f
(
t
)
=
(
C
υ
y
+
32
C
2
)
(
1
−
e
−
t
/
C
)
−
32
C
t
x
=
r
(
t
)
=
C
υ
x
(
1
−
e
−
t
/
C
)
\begin{array}{l} y = f(t) = (C{\upsilon _y} + 32{C^2})\left( {1 - {e^{ - t/C}}} \right) - 32Ct\\ x = r(t) = C{\upsilon _x}\left( {1 - {e^{ - t/C}}} \right) \end{array}
y=f(t)=(Cυy+32C2)(1−e−t/C)−32Ctx=r(t)=Cυx(1−e−t/C)
b
0
=
4
5
0
,
υ
y
=
υ
x
=
160
f
t
/
sec
,
C
=
10
{b_0} = {45^0},{\upsilon _y} = {\upsilon _x} = 160ft/\sec ,C = 10
b0=450,υy=υx=160ft/sec,C=10
f
(
8
)
=
83.22
,
f
(
9
)
=
−
31.5
f(8) = 83.22,f(9) = - 31.5
f(8)=83.22,f(9)=−31.5
p
0
=
8
,
f
′
(
p
0
)
=
−
104.32
{p_0} = 8,f'({p_0}) = - 104.32
p0=8,f′(p0)=−104.32
g
(
x
)
=
x
−
f
(
x
)
f
′
(
x
)
⇒
p
1
=
8
−
83.22
−
104.32
=
8.7977
g(x) = x - \frac{{f(x)}}{{f'(x)}} \Rightarrow {p_1} = 8 - \frac{{83.22}}{{ - 104.32}} = 8.7977
g(x)=x−f′(x)f(x)⇒p1=8−−104.3283.22=8.7977
可知炮弹落地时间位于8s~9s内，于是利用迭代的算法进行求解，其核心算法代码如下
void NewtonRaphsonIterEquation(float (*f)(float, float, float), float (*df)(float, float, float), float (*ddf)(float, float, float), float C, float t0, float Vy, float Vx, float** valset,int& nIter,float& yDis,float& xt, const int nMaxIter,const float eps,bool& ifOutPre)
{
nIter=0;
float pk=t0,pk1,err;
//solve out the fly time
do
{
pk1 = NR_Equation(C,Vy,pk, f, df);
valset[nIter][0] = nIter; valset[nIter][1] = pk; valset[nIter][2] = pk1-pk;
valset[nIter][3] = f(C,Vy,pk);
++nIter;
err = (pk1 - pk) / pk1;
pk = pk1;
if (nIter >= nMaxIter)
{
cout << "cannot get the root within the precision" << endl;
ifOutPre = true;
break;
}
} while (fabs(err) > eps);
//return the flytime
xt = pk1;
//solve out the y(t) derivative =0 time
int tIter=0;
pk = t0;
do
{
pk1 = NR_Equation(C, Vy, pk, df, ddf);
++tIter;
err = (pk1 - pk);
pk = pk1;
if (nIter >= nMaxIter)
{
break;
}
} while (fabs(err) > eps);
//the highest height
yDis = f(C, Vy, pk1);
}
5、实验数据结果
实验具体数据结果在文末链接中，大家可以自行下载验证。其中，各个算法的计算结果均正确，且算法收敛速度依次为 牛顿-拉弗森方法 > 试值法 > 二分法。因此可以看出，根据牛顿-拉弗森方法的原理，其收敛速度还与方程解附近的函数斜率有关，并不是所有情况下其收敛速度都会较快，具有严重的依赖性。当函数性质不太良好时，可能其收敛速度会低于两外两种方法。因此，在实际运用中，还应考虑函数的特性。
三、结论
通过分析求解非线性方程近似解的原理，我们发现可以用三种不同的迭代方式，使得方程的根的可能取值区间无限的缩小，最后可以通过比较取值区间的横向长度，或区间中点函数值与0的距离，如果在我们设置的允许误差范围内，并可认为找到了非线性方程的近似解。通过分析求解非线性方程近似解的原理，我们发现可以用三种不同的迭代方式，使得方程的根的可能取值区间无限的缩小，最后可以通过比较取值区间的横向长度，或区间中点函数值与0的距离，如果在我们设置的允许误差范围内，并可认为找到了非线性方程的近似解。
通过实验验证三种不同的方式，我们发现。三种方法求解非线性方程的根的收敛性质良好，且求得根的精确度非常高。
声明：本系列《数值计算方法》所有章节源代码均为作者编写，大家可以自行下载，仅供学习交流使用，欢迎大家评论留言或者私信交流，请勿用于任何商业行为。其中，各章实验数据结果也在资源链接中，大家可以自行实验对照，欢迎大家批评指正！文件中每个算法均由作者进行初步测试，部分算法可直接用于工程开发作为模块参考。
扫一扫
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
