《科学计算引论》实验课 第五章 非线性方程数值解法

1、科学计算引论科学计算引论实验实验第五章第五章 非线性方程数值解法非线性方程数值解法本章内容简介本章内容简介 5.1 几何方法几何方法 5.2 Picard 迭代法迭代法 5.3 Newton 迭代法迭代法F (X) = 0Matlab解方程函数解方程函数roots(p)：多项式的多项式的所有零点所有零点，p 是多项式系数向量是多项式系数向量fzero(f,x0)：求求 f(x)=0 在在 x0 附近的一个根，附近的一个根，f 是函是函数句柄，可以通过内联函数，匿名函数或函数文件来定义，数句柄，可以通过内联函数，匿名函数或函数文件来定义，但不能是方程或符号表达式！但不能是方程或符号表达式！sol

2、ve(f,v)：求方程关于指定自变量的解，求方程关于指定自变量的解，f 是是符号表达符号表达式式或或符号方程符号方程；l solve 也可解方程组也可解方程组 (包含非线性包含非线性)l 得不到解析解时，给出数值解得不到解析解时，给出数值解linsolve(A,b)：解线性方程组解线性方程组二分法二分法q 二分法的计算流程：二分法的计算流程：二分法二分法q 二分法的计算流程：二分法的计算流程：算法算法5.1 二分法二分法二分法示例二分法示例解解：弦截法、弦截法、Steffensen 方法方法q 弦截法：弦截法：111()()(, 1,2,.)kkkkkkkxxxxf xkf xf x 当当 时

3、，时，*kxx1()0kkkxxf x 需两个初始启动值需两个初始启动值q Steffensen方法：方法：改造弦截法得改造弦截法得Steffensen方法方法：21(), 1,2,.()()kkkkkkfxxxkf xf xf x Picard 迭代法迭代法q Picard迭代法迭代法基本思想基本思想l 构造构造 F (X) = 0 的一个等价方程组：的一个等价方程组： ()XX 得到一个迭代序列：得到一个迭代序列： X0，X1，X2，. . . 其中：其中：l 任取一个迭代初始值任取一个迭代初始值 X0 ，计算，计算 1()kkXX k = 0, 1, 2, . ( )( )( )12(,

4、)kkkTknXxxx 问题非线性方程非线性方程F (X) = 01212(,) , ()(),(),()TTnnnXxxxF XfXfXfXnXDl 迭代终止准则：迭代终止准则： （预设精度）（预设精度）1kkXX 误差估计误差估计定理定理5.1：设：设 于闭凸集于闭凸集 上连续上连续可可微微（G-可微），且满足可微），且满足这里这里 是从属于是从属于 中某向量范数中某向量范数 的矩阵范数，的矩阵范数， 则方程组则方程组 有有唯一不动点唯一不动点 ，迭代格式，迭代格式 自任意迭代初值自任意迭代初值 出发出发收敛收敛于该不动点，且有于该不动点，且有先验误差估计先验误差估计及及后验误差估计后验误

5、差估计000():sup()1, ,XDDDqX :nnDRR 0DD *101,2 ,1kkqXXXXqk *11,2, ,1kkkqXXXXqk nR()XX 00XD 1()kkXX *0XD 注：注： q (局部考察局部考察q*)越小，迭代收敛越快。越小，迭代收敛越快。加速迭代法加速迭代法Aitken 加速迭代公式加速迭代公式1()kkxx 11(),kkxx 2111111(),2kkkkkkkxxxxxxx q 方程组的加速迭代公式：方程组的加速迭代公式：11() ()kkkXIPXPX k = 0, 1, 2, . . q 标量方程的加速迭代公式：标量方程的加速迭代公式：Newt

6、on迭代法迭代法f (x) = 0111(), 1,2,.()kkkkf xxxkfx q Newton迭代法：迭代法：l 迭代终止准则：迭代终止准则： （预设精度）（预设精度） kkXX1| 或或 （预设精度）（预设精度） kF X|()| 11()(), 0,1,2,kkkkXXFXF Xk q 非线性方程组的非线性方程组的Newton迭代法迭代法：F (X) = 0简化简化Newton法法 & Newton下山法下山法q 简化的简化的Newton迭代法迭代法q Newton下山法下山法101(), 0,2)1, ,(kkkXXFFXXk l基本思想：减少导数的计算基本思想：减少导

7、数的计算 11()(),0,1,2,.kkkkXXFXF Xk l 下山因子的取法：下山因子的取法： 从从 =1 开始，逐次减半，直到满足下山条件开始，逐次减半，直到满足下山条件l 基本思想：要求每一步迭代满足基本思想：要求每一步迭代满足下山条件下山条件 1kkF XF X :01 l 具体做法：加具体做法：加下山因子下山因子保证全局收敛保证全局收敛简化简化Newton法法 & Newton下山法下山法q 简化的简化的Newton迭代法迭代法q Newton下山法下山法101(), 0,2)1, ,(kkkXXFFXXk l基本思想：减少导数的计算基本思想：减少导数的计算 11()()

8、,0,1,2,.kkkkXXFXF Xk l 下山因子的取法：下山因子的取法： 从从 =1 开始，逐次减半，直到满足下山条件开始，逐次减半，直到满足下山条件l 基本思想：要求每一步迭代满足基本思想：要求每一步迭代满足下山条件下山条件 1kkF XF X :01 l 具体做法：加具体做法：加下山因子下山因子保证全局收敛保证全局收敛Jacobi矩阵的计算矩阵的计算q 计算计算Jacobi矩阵矩阵的函数的函数% 求函数求函数fun在在x处的处的Jacobi矩阵矩阵function J=Jacobi(fun,x)m=length(x); % 确定变量的维度确定变量的维度n=length(fun(x);

9、 % 确定函数的维度确定函数的维度J=zeros(m,n); % 初始化初始化Jacobi矩阵矩阵epsnew=1.d-12; % 小增量小增量I=eye(m,n); % 保存单位向量的矩阵保存单位向量的矩阵for k=1:n % 依次计算依次计算Jacobi矩阵的第矩阵的第k列列 J(:,k)=(fun(x+epsnew*I(:,k)-fun(x)/epsnew;end例：例： f=(x) sin(x(1);cos(x(2); J=jacobi(f,0;0)标量方程的标量方程的Newton迭代算法迭代算法% fun和和dfun以字符串的形式输入即可以字符串的形式输入即可function xv

10、ect,xdif,fx,nit=newton(x0,nmax,tol,fun,dfun)err=tol+1; nit=0; xvect=x0; x=x0; fx=eval(fun); xdif=;while (nit tol), nit=nit+1; x=xvect(nit); dfx=eval(dfun); if (dfx = 0), nit=nmax; disp(Stop for vanishing dfun); else, xn=x-fx(nit)/dfx; err=abs(xn-x); xdif=xdif; err; x=xn; xvect=xvect;x; fx=fx;eval(fun); end;end;例：例： xvect,xdif,fx,ni