【非线性方程求解】二分法

1. 概述

介绍二分法求解含单变量非线性方程的过程。

2. 原理步骤

2.1 原理

现实中会遇到很多无法用线性等式描述的方程，以单变量为例，输入$x$后经过一系列复杂计算得到对应y，只能用$\phi \left(x\right)=y$表示映射关系，问题来了：当我们想要某个结果$y_m$ ，如何确定输入$x$，二分法可以帮你完成这一过程。
引入新的映射关系:$f\left(x\right)=\phi \left(x\right)-y_m$，表示的是随便给一个$x$，经过原先的复杂计算后，原输出与指定输出相差多少。这样原问题的求解就变成了求解非线性方程 $f\left(x\right)=0$。
求解非线性方程有很多方法，这里以最简单的二分法为例入手，通常实际中遇到的问题并不是单变量的，后续再来探讨非线性方程组的问题。

2.2 步骤

步骤如下：

1. 设置端点a,b和误差限 $\epsilon$；
2. 计算最大迭代次数k：
   $$k\ge \mathrm{ceil}\left(\frac{\ln (b-a)-\ln \epsilon }{\ln 2}\right)$$
   ceil表示向上取整；
3. 对$n=1,\dots ,k$循环执行d)和f)；
4. $c=(a+b)/2$，计算 $f\left(c\right)=0$
5. 若$f(c)\cdot f(b)<0$，则 $a=c$，否则$b=c$；
6. 当$|f(c)|<\epsilon$则停止循环，解为$c$。

3. 结论

介绍了二分法的原理和步骤。

4. 详细代码

测试

double MyFunction(double x) {
	return (x - 2)*(x - 2) - 4;
}
int main()
{
	CBinarySearch BinarySearch(100);
	try {
		BinarySearch.Search(-0.2, 1, 0.01, MyFunction);
		auto Record = BinarySearch.getRecord();
		cout << "i=" << Record.index << endl;
		cout << "x=" << Record.data[Record.index].c << endl;
		cout << "f=" << Record.data[Record.index].fc << endl;
	}
	catch (invalid_argument &e) {
		cout << e.what() << endl;
		//......
	}
	catch (out_of_range &e) {
		cout << e.what() << endl;
		//......
	}
	
	return 0;
}

CBinarySearch.h

/*
* Auther:	 sanfan66
*/
#pragma once
//#include <stdexcept>//exception

class CBinarySearch
{
public:
	struct SData {
		double a, b, c, fa, fb, fc;
	};
	struct SRecord {
		SData *data;
		int maxNum;
		int index;
	};

	CBinarySearch(int mMaxNum) {
		Record.maxNum = mMaxNum;
		Record.data = new SData[Record.maxNum]{};
	}
	~CBinarySearch() {
		delete[] Record.data;
		Record.data = nullptr;
	}
	double Search(double a, double b, double ObjFunEps, double(*ObjFunction)(double)) ;
	SRecord getRecord() {
		return Record;
	}

private:
	SRecord Record;
};

CBinarySearch.cpp

#include "CBinarySearch.h"
#include <iostream>//cout
#include <sstream>//stringstream
using namespace std;//cout
double CBinarySearch::Search(double a, double b, double ObjFunEps, double(*ObjFunction)(double)) throw(exception)
{
	//Record.maxNum = ceil((log(b-a)-log(ObjFunEps))/log(2));
	double c = (a + b) / 2;
	double fc = (*ObjFunction)(c);
	Record.data[0] = {a, b, c, (*ObjFunction)(a), (*ObjFunction)(b), fc};
	if (abs(Record.data[0].fc) < ObjFunEps) {
		Record.index = 0;
		return 0;
	}
	else if (Record.data[0].fa*Record.data[0].fb>0) {
		throw invalid_argument("二分法初值不合理");
	}
	Record.index = -1;//判断是否有解
	for (int i = 1; i< Record.maxNum; i++) {
		Record.data[i] = Record.data[i-1];

		if (Record.data[i].fc*Record.data[i].fb<0) {
			Record.data[i].a = Record.data[i].c;
			Record.data[i].fa = (*ObjFunction)(Record.data[i].a);
		}
		else {
			Record.data[i].b = Record.data[i].c;
			Record.data[i].fb = (*ObjFunction)(Record.data[i].b);
		}

		Record.data[i].c = (Record.data[i].a + Record.data[i].b) / 2;
		Record.data[i].fc = (*ObjFunction)(Record.data[i].c);

		if (abs(Record.data[i].fc) < ObjFunEps) {
			Record.index = i;
			break;
		}
	}

	if (Record.index<0) {
		stringstream throwString;
		throwString << "二分法迭代次数超过设定值Record.maxNum="<< Record.maxNum;
		throw out_of_range(throwString.str());
		//cout << "二分法经过"<< Record.maxNum <<"次迭代，未找到解" << endl;
		//abort();
	}

	return 0;
}