数值分析——二分法和牛顿迭代(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

智能推荐

python-简单用户交互小程序-程序员宅基地

文章浏览阅读872次。1、第一种方式:root@kali:~/python# vim userinput.py#!/usr/bin/python#--*-- coding:utf-8 --*--Name = raw_input('please your name:\n')Age = raw_input('please your age:\n')Sex = raw_input('please

在Redhat9 Linux下安装汉化eclipse3.1.2的c/c++开发平台-程序员宅基地

文章浏览阅读6.1k次。在Redhat9 Linux下安装汉化eclipse3.1.2的c/c++开发平台By: 吴垠Date: 2006-06-09Email: [email protected]版权信息:该文章版权由Wu Yin所有。可在非商业目的下任意传播和复制。对于商业 目的下对本文的任何行为需经作者同意。联系方式:lazy_fox#msn.com1. 本文需要的资源都可

atl dll中调用wtl_dll中使用wtl开发窗口-程序员宅基地

文章浏览阅读1.6k次。环境 xp sp3,vs20031.新建 ATL-ATL项目,项目名test0507,不要选中属性化,不要选中支持 mfc,2.stdafx.h添加#include #include #include #include #include #include #include #include 3.在原有的Ctest0507Module类基础上添加Add_dll中使用wtl开发窗口

【404 笔记】总结-程序员宅基地

文章浏览阅读174次。mvc mvp mvvmhttps://juejin.cn/post/6914971761163370509https://juejin.cn/post/6846687603547176974Handler谈谈消息机制Handler作用 ?有哪些要素 ?流程是怎样的 ?负责跨线程通信,这是因为在主线程不能做耗时操作,而子线程不能更新UI,所以当子线程中进行耗时操作后需要更新UI时,通过Handler将有关UI的操作切换到主线程中执行。https://www.jianshu.com/p/70d57

ZF网络架构深度详解-程序员宅基地

文章浏览阅读1w次,点赞50次,收藏121次。前言:ZF网络是2013年提出的,网上有很多关于它的介绍和讲解,但是很多内容讲的不太好(个人感觉),于是花时间收集了一些资料,整理了一些比较好的文章,从头到尾把ZFNet说了一遍。一、ZFNet简介1.1 为什么起名ZFnetwork ILSVRC 2013获胜者是来自Matthew Zeiler和Rob Fergus的卷积网络。它被称为ZFNet(Zeiler&..._zf网络

第二十章 注解-程序员宅基地

文章浏览阅读2.1w次,点赞29次,收藏44次。文章目录0.概述1.基本语法定义注解元注解2.编写注解处理器注解元素3.使用apt处理注解4.将观察者模式用于apt5.基于注解的单元测试让测试类继承自泛型类的一个特定版本0.概述注解(也被称为元数据)为我们在代码中添加信息提供了一种形式化的方法,使我们可以在稍后某个时刻非常方便地使用这些数据注解使我们能够以将由编译器来测试和验证的格式,存储有关程序的额外信息注解的优点:更加干净易读的代码以及编译器类型检查等@Override,表示当前的方法定义将覆盖超类中的方法。如果你不小心拼写错误,或者方

随便推点

visio将数据显示到形状上_visio形状数据-程序员宅基地

文章浏览阅读2.3k次。遇到的问题:“此数据图形无法应用到选定范围中的一个或多个无效形状”目的:将数据显示到形状上步骤: 1、定义形状的属性,即定义形状包含哪些数据;2、定义将哪些属性显示到形状上。将形状拖到画布中,右键,选择“数据”-&gt;定义形状数据。此处自定义形状具有的属性,如:设备名称、设备IP、设备位置等形状数据定义完成后,选择“编辑数据图形”,定义将哪些属性显示到形..._visio形状数据

STM32F4X ADC_stm32注入触发转换_hwx1546的博客-程序员宅基地

文章浏览阅读110次。ADC全称是Analog-Digital-Converter,模拟数字转换,也叫模数转换。为什么嵌入式系统需要ADC,我们知道在自然界中广泛存在模拟量,比如声、光、电、磁等。但是对于嵌入式系统来说,如果想要识别自然界的模拟量,就必须将模拟量转换成MCU熟悉的数字量,而在模拟量和数字量之间进行转换的模块就叫模数转换器。_stm32注入触发转换

侧滑菜单DrawerLayout+Toolbar结合使用,切换菜单-程序员宅基地

文章浏览阅读3.4k次,点赞2次,收藏6次。因为最近项目中用到侧滑菜单,于是就想到了谷歌提供的这个类网上看了看资料完成了侧滑菜单的功能,并写了一个demo记录下来,加深记忆和以后做类似功能时,直接看看笔记比住较方便主布局文件分为两个一个是侧滑菜单覆盖toolbar一个是不覆盖toolbar具体看效果图,无赖不会上传gif1.覆盖toolbar2、不覆盖toolbar2、menu主布局文件(

Boost Asio异步发送数据(async_write)崩溃问题记录_boost udp异步崩溃-程序员宅基地

文章浏览阅读3.3k次,点赞2次,收藏2次。背景服务端与客户端之间的网络通信(使用Boost Asio库异步编程模式实现),客户端会向服务端请求数据。在刚开始的测试中,是没有出现问题的。后来有一次测试时,服务端查询完数据后,向客户端发送时总是崩溃。通过gdb调试,可以发现是在调用到异步发送函数(boost::asio::async_write)后崩溃的。打印的栈信息如下:Program terminated with signa..._boost udp异步崩溃

CG动画制作项目第十一篇:后期剪辑以及音效、配音处理(一)-程序员宅基地

文章浏览阅读3.3k次。在做完所有的渲染过程之后,我们的工作相应的也进入到了后期阶段。后期需要处理的事情包括:视频剪辑、调色、特效处理、音效添加、配音、声音处理、BGM、二维镜头调整等等视频剪辑比较简单,目的是为了可以将一个完整流畅的动画展现在大家面前。可是我们如果单纯的拼接一组镜头,我们会感觉到这个动画显得非常的单调,然后我们就需要对我们的动画进行润色。其中之一是调色,调色对我们的动画的表现力有很大的提升,比如一个很悲...

帝国cms多表连接查询及结果集的调用-程序员宅基地

文章浏览阅读1.8k次。SQL语句:select //查询 a.*, //a表的所有列 b.truename as writer, //b表的truename 起个别名叫 writer ,实际对应a表中的字段w