目录

一.误差理论误差来源误差基本概念绝对误差相对误差有效数字误差传播应注意的问题MATLAB数值计算精度数据显示格式运算精度举例二. 非线性方程（组）求解二分法试位法简单迭代法牛顿法牛顿下山法割线法实例

一.误差理论

误差来源

1.模型误差

2.观测误差

3.截断误差：求解模型所用数值方法为近似方法因此有误差存在。

4.舍入误差

误差基本概念

绝对误差

给定一实数x∗x^*x∗,其近似值为xxx，则称e=x∗−xe=x^*-xe=x∗−x为近似值xxx的绝对误差,若存在一个正实数\varepsilon，使得∣e∣=∣x∗−x∣≤ε|e|=|x^*-x|\le\varepsilon∣e∣=∣x∗−x∣≤ε则称其为近似值的绝对误差限.

相对误差

相对误差公式：er=x∗−xx∗e_r=\frac{x^*-x}{x^*}er​=x∗x∗−x​

实用相对误差公式：erˉ=x∗−xx\bar{e_r}=\frac{x^*-x}{x}er​ˉ​=xx∗−x​

有效数字

误差传播

函数运算中误差的传播可以用函数的泰勒展开式进行估计：

f(x∗)−f(x)=f′(x)(x∗−x)+12f′′(ξ)(x∗−x)2f(x^*)-f(x)=f'(x)(x^*-x)+\frac{1}{2}f''(\xi){(x^*-x)}^2f(x∗)−f(x)=f′(x)(x∗−x)+21​f′′(ξ)(x∗−x)2

∣f(x∗)−f(x)∣≤∣f′(x)∣e(x)+12∣f′′(ξ)∣(e(x))2|f(x^*)-f(x)|\le|f'(x)|e(x)+\frac{1}{2}|f''(\xi)|{(e(x))}^2∣f(x∗)−f(x)∣≤∣f′(x)∣e(x)+21​∣f′′(ξ)∣(e(x))2

忽略e(x)e(x)e(x)的高阶项，其误差限为：

ε(f(x))≈∣f′(x)∣e(x)\varepsilon(f(x))\approx|f'(x)|e(x)ε(f(x))≈∣f′(x)∣e(x)

其解的相对误差是：

εrˉ(f(x))≈∣f′(x)∣e(x)f(x)\bar{\varepsilon_r}(f(x))\approx|f\prime(x)|\frac{e(x)}{f(x)}εr​ˉ​(f(x))≈∣f′(x)∣f(x)e(x)​

相对误差限为：

εr(f(x))=∣f′(x)∣e(x)f(x)\varepsilon_r(f(x))=|f\prime(x)|\frac{e(x)}{f(x)} εr​(f(x))=∣f′(x)∣f(x)e(x)​

应注意的问题

为了减小舍入误差的影响：

1.避免两个详尽的数相减

2.避免绝对值过小的数作除数

3.避免大数吃掉小数

4先化简再计算

5.采用数值稳定性好的算法

MATLAB数值计算精度

数据显示格式

控制数据输出格式的指令是format

format 数显标识符

运算精度

MATLAB提供了三种运算精度的算法：

数值算法：每个数值都取16位有效数字，运算速度最快。

符号算法：每个数据均为符号，精确但是慢，占用内存多。

可控精度算法：利用控制精度指令digits（n）使运算以n位有效数字进行，亦可用vpa（）函数

举例

秦九韶算法

function y=QJS(fun,x)fun=sym(fun);a=sym2poly(fun);s(1)=a(1);for k=2:length(a)s(k)=s(k-1)*x+a(k);endy=s(end);

二. 非线性方程（组）求解

二分法

优点：简单，对f（x）要求不高，连续即可

缺点：无法求复根和偶重根，收敛慢

function [x,fx,iter,X] =bisect(fun,a,b,eps,varargin)%二分法%二分法求非线性方程的根%输入参数：%——fun：待求根方程的函数%——a，b：初始区间的端点%——eps：精度要求，默认值为1e-6%——p1,p2,...：求根函数附加参数%输出参数：%——x:非线性方程组的近似根%——fx:根x处的函数值%——iter：迭代次数%——X：每一步迭代的结果if nargin<3error('输入参数至少需要3个!');endif nargin<4isepmty(eps);eps=1e-6;endfa=feval(fun,a,varargin{:});fb=feval(fun,b,varargin{:});k=1;if fa*fb>0warning(['区间[',num2str(a),',',num2str(b),']可能没有根'])elseif fa==0x=a;fx=fa;elseif fb==0x=b;fx=fb;elsewhile abs(b-a)>epsx=(a+b)/2;fx=feval(fun,a,varargin{:});if fa*fx>0a=x;fa=fx;elseif fb*fx>0b=x;fb=fx;elsebreakendX(k)=x;k=k+1;end enditer=k;

试位法

function [x,fx,iter,X] =test_bit(fun,a,b,eps,varargin)%试位法求非线性方程的根%输入参数：%——fun：待求根方程的函数%——a，b：初始区间的端点%——eps：精度要求，默认值为1e-6%——p1,p2,...：求根函数附加参数%输出参数：%——x:非线性方程组的近似根%——fx:根x处的函数值%——iter：迭代次数%——X：每一步迭代的结果if nargin<3error('输入参数至少需要3个!');endif nargin<4isepmty(eps);eps=1e-6;endfa=feval(fun,a,varargin{:});fb=feval(fun,b,varargin{:});k=1;if fa*fb>0warning(['区间[',num2str(a),',',num2str(b),']可能没有根'])elseif fa==0x=a;fx=fa;elseif fb==0x=b;fx=fb;elsewhile abs(b-a)>eps%控制试位法结束条件x=b-(b-a)/(fa-fb)*fb;fx=feval(fun,x,varargin{:});%计算x处的函数值if fa*fx>0%条件a=x;%端点更新fa=fx;%端点函数值更新elseif fb*fx>0%条件b=x;%端点更新fb=fx;%端点函数值更新elsebreakendX(k)=x;k=k+1;end enditer=k;

简单迭代法

function [x,iter,X] = fixpt(g,x0,eps,maxiter)%简单迭代法求非线性方程的根%输入参数：%——g：迭代函数%——x0：初始迭代点%——eps：精度要求，默认值为1e-6%——maxiter：最大迭代次数，默认1e4%输出参数：%——x:非线性方程组的近似根%——iter：迭代次数%——X：每一步迭代的结果if nargin<2error('输入参数至少需要3个');endif nargin<3isempty(eps);eps=1e-6;endif nargin<4isempty(maxiter);maxiter=1e4;enderr=1;k=0;while err>epsk=k+1;x1=feval(g,x0);err=abs(x1-x0);x0=x1;X(k)=x1;endif k==maxitererror('迭代次数超限')endx=x1;iter=k;X=X';end

收敛才可以迭代出结果

牛顿法

function [x,iter,X] = newton(fun,x0,eps,maxiter)%牛顿法求非线性方程的根%输入参数：%——fun：迭代函数%——x0：初始迭代点%——eps：精度要求，默认值为1e-6%——maxiter：最大迭代次数，默认1e4%输出参数：%——x:非线性方程组的近似根%——iter：迭代次数%——X：每一步迭代的结果if nargin<2error('输入参数至少需要3个');endif nargin<3isempty(eps);eps=1e-6;endif nargin<4isempty(maxiter);maxiter=1e4;enderr=1;k=0;while err>epsk=k+1;[fx0,dfx0]=feval(fun,x0);if dfx0==0error('f（x）在x0处的导数为0，停止计算')endx1=x0-fx0/dfx0;err=abs(x1-x0);x0=x1;X(k)=x1;endif k>=maxitererror('迭代次数超限')endx=x1;iter=k;X=X';end

牛顿下山法

function [x,iter,X] = newton_down(fun,x0,eps,maxiter)%牛顿下山法求非线性方程的根%输入参数：%——fun：迭代函数%——x0：初始迭代点%——eps：精度要求，默认值为1e-6%——maxiter：最大迭代次数，默认1e4%输出参数：%——x:非线性方程组的近似根%——iter：迭代次数%——X：每一步迭代的结果if nargin<2error('输入参数至少需要3个');endif nargin<3isempty(eps);eps=1e-6;endif nargin<4isempty(maxiter);maxiter=1e4;end[fx0,dfx0]=feval(fun,x0);k=0;tol=fx0/dfx0;while abs(tol)>epslambda=1;[fx0,dfx0]=feval(fun,x0);x1=x0-lambda*fx0/dfx0;[fx1,dfx1]=feval(fun,x1);while abs(fx1)>=abs(fx0)%保证迭代后的函数值比之前小lambda=lambda/2;x1=x0-lambda*fx0/dfx0;%牛顿迭代下山[fx1,dfx1]=feval(fun,x1);x0=x1;endtol=fx1/dfx1;k=k+1;x0=x1;X(k)=x1;endif k>=maxitererror('迭代次数超限')endx=x1;iter=k;X=X';end

割线法

function [x,iter,X] = secant(fun,x0,eps,maxiter)%割线法求非线性方程的根%输入参数：%——fun：迭代函数%——x0：初始迭代点%——eps：精度要求，默认值为1e-6%——maxiter：最大迭代次数，默认1e4%输出参数：%——x:非线性方程组的近似根%——iter：迭代次数%——X：每一步迭代的结果if nargin<2error('输入参数至少需要3个');endif nargin<3isempty(eps);eps=1e-6;endif nargin<4isempty(maxiter);maxiter=1e4;endf0=feval(fun,x0);%计算x0处的函数值f1=feval(fun,x1);%计算x1处的函数值k=0;while err>epsk=k+1;x2=x1-f1*(x1-x0)/(f1-f0);f0=f1;x0=x1;x1=x2;f1=feval(fun,x1);%计算x1处的函数值X(k)=x1;endif k>=maxitererror('迭代次数超限')endx=x2;iter=k;X=X';end

实例

求解humps函数在区间[0,2]的零点

x=0:0.01:2;fun=@(x)1./((x-.3).^2+.01)+1./((x-.9).^2+.04)-6;%定义函数表达式[x1,fx1]=bisect(fun,0,2,1e-10);%二分法求解plot(x,humps(x),[0,2],[0,0],'k',x1,fx1,'k*')text(1,40,['humps函数零点{\itx}^*=',num2str(x1)])title('humps函数零点求解')