文章目录
- 图解法
- 一般步骤
- 实现代码
- 实例分析
- 二分法
- 一般步骤
- 实现代码
- 实例分析
- 试位法
- 实现代码
- 实例分析
求解非线性方程的根是科学研究和工程计算中另一个非常重要的问题. 无论是在数学理论上的微分方程求解, 优化问题, 数值模拟, 还是在工程科学中的诸如流体动力学, 结构力学, 电磁场等领域, 以及在一些经济学科和其他社会科学学科中的模型研究中, 经常会遇到求解非线性方程.
非线性方程的数值解法大致分为划界法与迭代法两类. 每一类方法中都有很多经典的方法, 各有特色. 划界法主要有图解法, 二分法以及试位法; 而迭代法主要有不动点迭代法, 牛顿法, 割线法, 抛物线法以及逆二次插值法. 此外, 还有结合了两者的布伦特法.
本文代码所需的预处理如下:
#include <QStandardPaths>
#include <functional>
#include <QFile>
#include <math.h>
using namespace std;
由于本文代码是在Qt项目中写的, 使用了QFile
, 读者也可以将其改为使用标准库的文件流实现.
图解法
一般步骤
- 分别绘制曲线 y = f ( x ) y=f(x) y=f(x)和直线 y = 0 y=0 y=0的图形.
- 观察 y = f ( x ) y=f(x) y=f(x)和 y = 0 y=0 y=0的交点, 该交点的横坐标即为方程的根, 该值只能在图上读出(必要时可利用放大工具将图形局部放大), 因此一般情况下精度不会很高.
实现代码
/*
* 导出画图数据(一元函数)
* f : 函数
* a : 画图区间左端点
* b : 画图区间右端点
* n : 区间划分数
* path : 文件保存路径
* title : 表头
*
* 返回(bool):
* true :导出失败
* false :导出成功
*/
bool plot(function<double(double)> f, double a = 0, double b = 1, unsigned n = 1000, QString path = QStandardPaths::writableLocation(QStandardPaths::DesktopLocation), const char *title = "x,f(x)")
{
QFile file(path);
if (!file.open(QIODevice::WriteOnly))
return true;
if (b < a)
{
double t(a);
a = b;
b = t;
}
file.write(title);
double step((b - a) / n++);
do
{
double f_a(f(a));
if (isnan(f_a))
{
file.close();
return true;
}
file.write(('\n' + QString::number(a) + ',' + QString::number(f_a)).toUtf8());
a += step;
} while (--n);
return false;
}
实例分析
求函数 f ( x ) = e x + x − 2 f(x)=e^x+x-2 f(x)=ex+x−2在 [ 0 , 2 ] [0,2] [0,2]上的根.
绘制函数图像如下:
从图中可以看出, 方程的根约为
0.45
0.45
0.45, 而利用MATLAB求得的精确解为
0.4428544010
⋯
0.4428544010\cdots
0.4428544010⋯, 误差在允许的范围内, 但若要寻求更高的精度则需要其它数值方法.
二分法
二分法是一种用于求解非线性方程的数值方法. 它的基本思想是将函数的定义域分为两个区间, 然后在每个区间内使用线性近似来寻找方程的根. 通过不断缩小搜索区间, 直到找到满足精度要求的解.
在二分法中, 我们首先确定一个初始区间, 该区间包含方程的解. 然后, 我们计算该区间的中点, 并计算中点处的方程值. 如果中点处的方程值小于零, 则解位于中点右侧的子区间内; 如果中点处的方程值大于零, 则解位于中点左侧的子区间内. 我们将根据这个信息, 将原始区间划分为两个新的子区间, 其中一个包含解, 另一个不包含解.
重复上述步骤, 每次将包含解的子区间缩小一半, 直到包含解的子区间足够小, 使得我们能够确定解的近似值.
一般步骤
输入: 实数区间 [ a , b ] [a,b] [a,b], 函数 f ( x ) f(x) f(x), 计算精度 ε \varepsilon ε
输出: 近似解 x x x
- 如果
f
(
a
)
f
(
b
)
>
0
f(a)f(b)>0
f(a)f(b)>0 则
- 返回 无解
- 循环,直到
b
−
a
<
ε
b-a<\varepsilon
b−a<ε 为止
- m ← a + b 2 m\gets\dfrac{a+b}2 m←2a+b
- 如果
f
(
m
)
f
(
a
)
>
0
f(m)f(a)>0
f(m)f(a)>0 则
- a ← m a\gets m a←m
- 否则
- b ← m b\gets m b←m
实现代码
/*
* 二分法
* f : 函数
* a : 区间左端点
* b : 区间右端点
* eps: 终止条件
*
* 返回(double): 近似解
*/
double dichotomy(const function<double(double)> &f, double a, double b, double eps)
{
if (a > b)
{
double t = a;
a = b;
b = t;
}
if (f(a) > 0)
{
if (f(b) > 0)
return NAN;
while (b - a > eps)
{
double m((a + b) / 2);
if (isnan(f(m)))
return NAN;
if (f(m) > 0)
if (a == m)
break;
else
a = m;
else
if (b == m)
break;
else
b = m;
}
}
else
{
if (f(b) < 0)
return NAN;
while (b - a > eps)
{
double m((a + b) / 2);
if (isnan(f(m)))
return NAN;
if (f(m) > 0)
if (b == m)
break;
else
b = m;
else
if (a == m)
break;
else
a = m;
}
}
return (a + b) / 2;
}
实例分析
利用二分法求函数 f ( x ) = sin x − x − x 2 + 0.5 f(x)=\sin x-x-x^2+0.5 f(x)=sinx−x−x2+0.5在区间 [ 0 , 1 ] [0,1] [0,1]内的零点.
代入程序可得区间左右端点变化情况如图所示:
试位法
试位法是二分法的一个变形. 在二分法中, 若把取中点改为取点
(
a
,
f
(
a
)
)
(a,f(a))
(a,f(a))与
(
b
,
f
(
b
)
)
(b,f(b))
(b,f(b))同
x
x
x轴的线性插值, 即
m
=
b
−
b
−
a
f
(
b
)
−
f
(
a
)
f
(
b
)
m=b-\dfrac{b-a}{f(b)-f(a)}f(b)
m=b−f(b)−f(a)b−af(b)
就得到了试位法.
实现代码
/*
* 二分法
* f : 函数
* a : 区间左端点
* b : 区间右端点
* eps: 终止条件
*
* 返回(double): 近似解
*/
double dichotomy(const function<double(double)> &f, double a, double b, double eps)
{
if (a > b)
{
double t = a;
a = b;
b = t;
}
if (f(a) > 0)
{
if (f(b) > 0)
return NAN;
while (b - a > eps)
{
double m(b - (b - a) * f(b) / (f(b) - f(a)));
if (isnan(f(m)))
return NAN;
if (f(m) > 0)
if (a == m)
break;
else
a = m;
else
if (b == m)
break;
else
b = m;
}
}
else
{
if (f(b) < 0)
return NAN;
while (b - a > eps)
{
double m(b - (b - a) * f(b) / (f(b) - f(a)));
if (isnan(f(m)))
return NAN;
if (f(m) > 0)
if (b == m)
break;
else
b = m;
else
if (a == m)
break;
else
a = m;
}
}
return b - (b - a) * f(b) / (f(b) - f(a));
}
实例分析
用试位法求解上例方程.
可得区间左右端点变化情况如图所示: