数值分析学习——从抄书开始

最新推荐文章于 2024-09-02 22:08:11 发布

我是西红柿呀

最新推荐文章于 2024-09-02 22:08:11 发布

阅读量1.4k

点赞数 5

分类专栏：数值分析

本文链接：https://blog.csdn.net/qq_40248871/article/details/89673704

版权

本文通过学习《现代数值分析(MATLAB版)》一书，探讨非线性方程的求根方法，包括二分法、简单迭代法、Aitken加速、牛顿法和割线法，并使用MATLAB、Python及R语言实现这些方法。旨在结合理论与实践，加深对数值分析的理解。

摘要由CSDN通过智能技术生成

前言

既然书中有了详尽的理论分析、公式推导以及源程序示例，还要抄书呢? 一方面督促自己学习, 另一方面也便于自己熟悉一些语言工具以及数值分析中的基本概念.

采用哪些工具

本次学习数值分析将会采用 MATLAB、Python 以及 R语言作为工具.

使用的教材

现代数值分析(MATLAB版) 马昌凤编著.
下面愉快地进入到第二章的学习吧!

非线性方程的求根方法

本章涉及到的求解方法本质上都是迭代法, 只是迭代函数在形式上的不同!

二分法

二分法的 MATLAB、Python、R程序分别如下

function [x,k] = mbisec(f,a,b,ep)
	%{
	Function: 二分法求 f(x)=0 在区间 [a,b] 中的一个根
	Input: 
		f: 输入函数
		a: 左端点
		b: 右端点
		ep: 精度
	Process: 二分法
	Output: 
		x: 近似解
		k: 迭代数
	%}
	x = (a + b) / 2.0; k = 0;
	while abs(feval(f,x)) > ep | (b - a) > ep
		if feval(f,x) * feval(f,a) < 0
			b = x;
		else
			a = x;
		end
		x = (a + b) / 2.0; k = k+1;
	end
end

%{
Example:
f = @(x) x - exp(-x);
[x,k] = mbisec(f,0,1,1e-5)
%}

import numpy as np


def mbisec(f,a,b,ep=1e-5):
	"""二分法求解函数的零点""" 
	#利用 Python 写的程序运行效率要比 MATLAB 好，打开MATLAB太慢了！
	x = (a+b) / 2
	k = 0
	while np.abs(f(x)) > ep or b-a > ep:
		if f(x) * f(a) < 0:
			b = x
		else:
			a = x
		x = (a + b) / 2
		k += 1
	print(x,k)

def my_f(x):
	return x - np.exp(-x)

#mbisec(f=my_f, a=0, b=1)

#page31T2.1
g = lambda x: np.exp(x) + x -4
mbisec(f=g, a=1, b=2, ep=1e-4)

mbisec <- function(f,a,b,ep=1e-5){
  #二分法求根
  x <- (a+b)/2
  k <- 0
  while(abs(f(x)) > ep || b-a >ep){
    if(f(x) * f(a) <0){
      b = x
    } else {
      a = x
    }
    x = (a+b)/2
    k <- k + 1
  }
  out=c(x,k)
  return(out)
}

my_fun <- function(x){
  x - exp(-x)
}

out = mbisec(f=my_fun,a=0,b=1)
print(out)

从语法上看, 三者没啥太大区别, 唯一不同的就是 MATLAB 的默认值设置有点麻烦.

简单迭代法

function [x, k] = miter(phi, x0, ep, N)
	%Function: 利用简单迭代法求解函数根
	%Input: 
	%    phi(x): 迭代函数
	%    x0: 初始点
	%    ep: 精度
	%    N: 最大迭代次数
	%Output:
	%   x: 所求根
	%   k: 迭代次数

	%默认值设置
	if nargin < 4 N = 500; end
	if nargin < 3 ep = 1e-4; end

	%初始化迭代次数
	k = 0;
	while k < N
		x = feval(phi, x0);
		if abs(x-x0) < ep break; end
		x0 = x;