数值分析——二分法和牛顿迭代(Bisection Method & Newton‘s Method)_二分法公式_怀帝阍而不见的博客-程序员秘密

技术标签: c++  计算数学  

本系列整理自博主21年秋季学期本科课程 数值分析I 的编程作业,内容相对基础,参考书: David Kincaid, Ward Cheney - Numerical Analysis Mathematics of Scientific Computing (2002, Americal Mathematical Society)

目录

适用对象

二分法

(1)算法描述

(2)代码

牛顿迭代法

(1)算法描述

 (2)代码

(3)定理

1. Theorem on Newton's Method

2. Theorem on Newton's Method for Convex Function

拓展

(1)弦方法

(2)函数迭代(Functional Iteration)

1. 压缩映射定理(Contractive Mapping Theorem)

2. 误差分析

总结


适用对象

求非线性方程的数值解,即寻找x使得目标函数f(x)=0

二分法

(1)算法描述

要求f(x)在区间[a,b]上连续,且f(a)f(b) < 0.

令c=a+(b-a)/2,检验f(a)f(c) < 0是否为真,若真,则令b=c;若否,且f(a)f(c)!=0,则令a=c.

重复上述过程,直到f(a)f(c)=0即f(c)=0时停止。数值上,一般要求|f(c)|的值小于某个数,或区间长度小于某个数时,停止程序。

该算法线性收敛于方程的根。

(2)代码

#include<iostream>
#include<math.h>
#include <iomanip>
using namespace std;
int M 50;//set the step M large enough
double f(double x){}//the target function

int sign(double x) {
	if (x > 0)
		return 1;
	else if (x == 0)
		return 0;
	else return -1;
}
int bisection(double a, double b) {
	//apply bisection on function f on interval [a,b]
	double left = a, right = b, length = mid = 0, valuel = f(left), valuer = f(right), valuem = 0;
	cout << "---------------------------------------------------------------------------" << endl;
	cout << setw(2) << " k" 
		 << setw(10) << "left"<< setw(13) << "f(left)" 
		 << setw(10) << "right" << setw(15) << "f(right)" 
		 << setw(10) << "mid" << setw(13) << "f(mid)" << endl;
	cout << "---------------------------------------------------------------------------" << endl;
	if (sign(valuel) == sign(valuer)) {
		length = right - left;
		mid = left + length / 2;
		valuem = f(mid); valuel = f(left); valuer = f(right);
		cout << setw(6) << 1 
			 << setw(10) << left<< setw(15) << valuel 
		     << setw(10) << right << setw(15) << valuer 
			 << setw(10) << mid << setw(15) << valuem<< endl;
		cout << "---------------------------------------------------------------------------" << endl;
		cout << "f(x) signs same at both ends of the interval ["< a < <", "< < b < <"], cannot use bisection algrithm to find the root." << endl;
		return -1;
	}
	for (int k = 1; k <= M; k++) {
		length = right - left;
		mid = left + length / 2;
		valuem = f(mid); valuel = f(left); valuer = f(right);
		cout << setw(6) << k
			 << setw(10) << left << setw(15) << valuel
			 << setw(10) << right << setw(15) << valuer
			 << setw(10) << mid << setw(15) << valuem << endl;
		if (abs(length) < 0.0001 || abs(valuem) < 0.0001)//stop when the value of function is smaller than 0.0001
			break;
		if (sign(valuem) != sign(valuel))
			right = mid;
		else
			left = mid;
	}
	cout << "---------------------------------------------------------------------------" << endl;
	if (abs(length) < 0.0001 && abs(valuem) > 100)//this determination condition should be modified according to different functions
		cout << "f(x) has no root on interval [" << a << ", " << b << "],and it tend to infinity at the point x=" << mid  << endl;
	else
		cout << "the root of f(x) on interval [" << a << "," << b << "] is x = " << mid << endl;
	return 1;
}

牛顿迭代法

(1)算法描述

令r为f的零点,x趋向于r,由Taylor定理:0=f(r)=f(x+h)=f(x)+h*f'(x)+O(h^2),其中O(h^2)可以忽略,推出h= - f(x)/f'(x),从而零点 r=x+h=x - f(x)/f'(x)

牛顿迭代公式:

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

该算法二次收敛于方程的根。

 (2)代码

#include<iostream>
#include<math.h>
#include <iomanip>
using namespace std;
int M = 50;//the maximum of operation step
double f(double x){}//the target function
double df(double x){}//the derivative of f

double Newton_Method(double x) {//x is the initial value
	cout << " Newton's method" << endl;
	cout << "-----------------------------------------------" << endl;
	cout << setw(9) << " k" << setw(13) << "x_k" 
         << setw(15) << "f(x_k)" << setw(15) << "df(x_k)"
		 << resetiosflags(ios::left) << endl;
	cout << "-----------------------------------------------" << endl;
    double y = f(x);
	for (int k = 1; k <= M; k++) {
		cout << setw(7) << k << setw(13) << x 
             << setw(15) << y << setw(15) << df(x)<< endl;
		if (abs(y) == 0)
			break;
		x = x - y / df(x);
		y = f(x);
	}
	cout << "-----------------------------------------------" << endl;
	cout << "the root of the function f is x = " << x << endl;
	return 0;
}

(3)定理

1. Theorem on Newton's Method

对于牛顿迭代,若 f'' 连续,r 为 f 单根,则存在 r 的领域 B 和常数 c,当初值 x0∈B 时,xn 将稳定趋于 r,且二阶收敛

\small |x_{n+1}-r|\leqslant c(x_n-r)^2

若 r 为 k 重根,将迭代公式修改为

\small x_{n+1}=x_n-\frac{kf(x_n)}{f'(x_n)}

且 xn 一阶收敛于 r

2. Theorem on Newton's Method for Convex Function

若 f 二次连续可导,单调,凸,且存在唯一零点,则任意初值,牛顿迭代都将收敛于该零点。

拓展

(1)弦方法

x_{n+1}=x_n-f(x_n)\frac{x_n-x_{n-1}}{f(x_n)-f(x_{n-1})}

超线性收敛:

r 为零点,记 en=xn-r,由Taylor定理:f(xn)=f(r+en)=f(r)+en f'(r)+1/2 en^2 f''(r) +O(en^3),推出 

f(xn)/en=f'(r)+1/2 en f"(r)+O(en^2),同理f(xn-1)/en-1=f'(r)+1/2 en-1 f"(r)+O(en-1^2)

由(xn-xn-1)/(f(xn)-f(xn-1))≈1/f'(r),

\small e_{n+1}=\frac{x_n-x_{n-1}}{f(x_n)-f(x_{n-1})}\frac{f(x_n)/e_n-f(x_{n-1})/e_{n-1}}{x_n-x_{n-1}}e_ne_{n-1} \approx \frac{1}{2}\frac{f''(r)}{f'(r)}e_ne_{n-1}

记C=1/2 f''(r)/f'(r),则  \small e_{n+1}=Ce_n e_{n-1}

注:

弦方法收敛阶高于二分法,低于牛顿法。相比于牛顿迭代,弦方法不需要对函数求导,用“弦”来代替导数。

(2)函数迭代(Functional Iteration)

x_{n+1}=F(x_n)

1. 压缩映射定理(Contractive Mapping Theorem)

设C为实数集R的闭子集,F:C->C压缩映射,则F有唯一不动点s,即任意初值x0∈C,xn收敛于s,且{xn}为Cauchy列。

2. 误差分析

记en=xn-s,若F'存在且连续,由中值定理\small x_{n+1}-s=F(x_n)-F(s)=F'(\zeta _n)(x_n-s),即

\small e_{n+1}=F'(\zeta _n)e_n

在不动点s处,设\small F^{(k)}(s)=0,1\leqslant k<q,且\small F^{(q)}(s)\neq 0,则\small e_{n+1}=F(x_n)-F(s)=F(s+e_n)-F(s)

对F(s+en)做Taylor展开,可得

\small e_{n+1}=\frac{1}{q!}F^{(q)}(\zeta_n)e_n^q

即迭代法 q阶收敛于s

总结

算法不足:由于函数f使用全局定义,故程序一次只能操作一个函数,不能同时对多个函数求根。

版权声明:本文为博主原创文章,遵循 CC 4.0 BY-SA 版权协议,转载请附上原文出处链接和本声明。
本文链接:https://blog.csdn.net/lk32767/article/details/122673034

智能推荐

大数定律与中心极限定理_白水baishui的博客-程序员秘密

大数定律定义:设X1,X2,...,Xn,...X_1,X_2,...,X_n,...X1​,X2​,...,Xn​,...为随机变量序列,XXX为随机变量,若对任意的正数ϵ\epsilonϵ有:lim⁡n→∞P(∣Xn−X∣⩾ϵ)=0\lim_{n\to\infty}P(|X_n-X|\geqslant\epsilon)=0n→∞lim​P(∣Xn​−X∣⩾ϵ)=0或lim⁡n→∞P(∣Xn...

[NEFU ACM大一暑假集训 解题报告]前缀和与差分_鱼竿钓鱼干的博客-程序员秘密

[NEFU ACM大一暑假集训 解题报告]前缀和与差分题量略大,所以解题报告和fjy大佬分了一下工由我负责A-K部分题解(不是AK部分题解啊,哈哈)后半部分题解(LM+R~V+XYZ)由fjy大佬发布剩下的题目烦请有余力的同学自行解决题目A - Molly’s Chemicals题目求区间和为m的幂次的子区间个数,也就是要满足下面这个公式。sum[r]−s[l−1]=mk(k&gt;=0)sum[r]-s[l-1]=m^k(k&gt;=0)sum[r]−s[l−1]=mk(k&gt;=0)

sourceinsight使用技巧_Jasonfang0118的博客-程序员秘密

1 sourceinsight screen font 的默认字体是Verdana的,它是一直变宽字体。在Document style中可以将字体改为定宽的Courier 2 勾掉indent Open Brace和Indent Close Brace的效果: 继上一段,在相对缩进行里, 如果输入"{"或"}", 则自动和上一行列对齐 3 今天把一个用sourceinsight排版整齐

ETL工具Kettle,下载与安装部署_toxcode的博客-程序员秘密

来源:微信公众号 -DD程序鹅原文:https://mp.weixin.qq.com/s/VKUy9mvzv28gNwPlU1X1qA版权声明:本文为博主原创文章,转载请附上原文链接!更多系列可以搜索上面公众号,提前查阅。上篇介绍了Kettle是什么、概念模型和核心组件,相信大家已经对Kettle有了初步认识。该篇主要介绍Kettle下载与window中安装部署。下载官方下载地址:https://community.hitachivantara.com/s/art...

Apache Spark 2.2.0 中文文档 - Spark 编程指南 | ApacheCN_spark2.2.0中文文档_片刻小哥哥的博客-程序员秘密

Spark 编程指南概述Spark 依赖初始化 Spark使用 Shell弹性分布式数据集 (RDDs)并行集合外部 Datasets(数据集)RDD 操作基础传递 Functions(函数)给 Spark理解闭包示例Local(本地)vs. cluster(集群)模式打印 RDD 的 elemen

分页逻辑实现_思维态度行动的博客-程序员秘密

分页简介用来将数据分割成多个部分来分页面展示。什么时候用?数据量达到一定的时候,就需要使用分页来进行数据分割。要不然可能会面临以下问题:客户端一次性显示太多数据会影响到用户的体验,比如很难找到客户想要的信息,以及加载页面数据过慢。对于服务端来说,一次性传送的数据过多,可能造成内存溢出。分页的分类分页分为真分页和假分页:真分页(物理分页):在mysql中使用select * ...

随便推点

MYSQL建表规则_奋斗的小虾米的博客-程序员秘密

建立表规约【强制】表名、字段名必须使用小写字母或数字,禁止出现数字开头,禁止两个下划线中间只 出现数字。数据库字段名的修改代价很大,因为无法进行预发布,所以字段名称需要慎重考虑。说明:MySQL 在 Windows 下不区分大小写,但在 Linux 下默认是区分大小写。因此,数据库名、表名、字段名,都不允许出现任何大写字母,避免节外生枝。正例:aliyun_admin,rdc_confi...

C# 常量分类与测试_Allenxixi的博客-程序员秘密

C# 常量分类与测试关注我,带你了解C#的魅力

FFmpeg API 之 libavcodec库_H&A的博客-程序员秘密

libavformat 库负责封装和解封装,而 libavcodec 则用于解码和编码。类型 AVPacket 表示编码后的数据,其中包含一个或多个编码后的帧数据。类型 AVFrame 表示解码后,或者说原始的帧数据。编码和解码在某种程度来说,就是两者之间的互相转换。1、编解码概述FFmpeg 提供的 encode/decode API 有如下四个函数 avcodec_send_packet () / avcodec_receive_frame () / avcodec_send_fram.

vue mescroll使用上拉加载、下拉刷新_黎轩栀海的博客-程序员秘密

一个页面多次调用mescroll注:this.navData[1].mescroll.resetUpScroll(); // 点击筛选重置方法&lt;template&gt; &lt;div class="healthy-bean-record-c"&gt; &lt;div class="top"&gt; [外链图片转存失败(img-pPYyrr3K-1562119870434)(https://mp.csdn.net/assets/images/mine/tixianji

freenas 体验_weixin_33698823的博客-程序员秘密

简介FreeNAS是一套免费的NAS服务器,它能将一部普通PC变成网络存储服务器。该软件基于FreeBSD, Python,支持CIFS (samba), FTP, NFS protocols, Software RAID (0,1,5) 及 web 界面的设定工具安装安装过程很简单,几乎都是下一步下一步,但是注意的是fre...

oracle调用plsql函数,oracle plsql开发之一:函数_weixin_39765290的博客-程序员秘密

ORACLE 提供可以把PL/SQL 程序存储在数据库中,并可以在任何地方来运行它。这样就叫存储过程或函数。过程和函数统称为PL/SQL子程序,他们是被命名的PL/SQL块,均存储在数据库中,并通过输入、输出参数或输入/输出参数与其调用者交换信息。过程和函数的唯一区别是函数总向调用者返回数据,而过程则不返回数据。创建函数的语法如下:Sql代码 CREATE[ORREPLACE]FUNCTIONfu...

推荐文章

热门文章

相关标签