科學(xué)計(jì)算引論實(shí)驗(yàn)課

1、科學(xué)計(jì)算引論實(shí)驗(yàn)第五章 非線(xiàn)性方程數(shù)值解法2本章內(nèi)容簡(jiǎn)介 5.1 幾何方法 5.2 Picard 迭代法 5.3 Newton 迭代法F (X) = 03Matlab解方程函數(shù)roots(p)：多項(xiàng)式的所有零點(diǎn)，p 是多項(xiàng)式系數(shù)向量fzero(f,x0)：求 f(x)=0 在 x0 附近的一個(gè)根，f 是函數(shù)句柄，可以通過(guò)內(nèi)聯(lián)函數(shù)，匿名函數(shù)或函數(shù)文件來(lái)定義，但不能是方程或符號(hào)表達(dá)式！solve(f,v)：求方程關(guān)于指定自變量的解，f 是符號(hào)表達(dá)式或符號(hào)方程； solve 也可解方程組 (包含非線(xiàn)性) 得不到解析解時(shí)，給出數(shù)值解linsolve(A,b)：解線(xiàn)性方程組4二分法 二分法的計(jì)算流程：5

2、二分法 二分法的計(jì)算流程：6算法5.1 二分法7二分法示例解：8弦截法、Steffensen 方法 弦截法：當(dāng) 時(shí)，需兩個(gè)初始啟動(dòng)值 Steffensen方法：改造弦截法得Steffensen方法：9Picard 迭代法 Picard迭代法基本思想 構(gòu)造 F (X) = 0 的一個(gè)等價(jià)方程組：得到一個(gè)迭代序列： X0，X1，X2，. . . 其中： 任取一個(gè)迭代初始值 X0 ，計(jì)算 k = 0, 1, 2, . 問(wèn)題非線(xiàn)性方程F (X) = 0 迭代終止準(zhǔn)則： （預(yù)設(shè)精度）10誤差估計(jì)定理5.1：設(shè) 于閉凸集 上連續(xù)可微（G-可微），且滿(mǎn)足這里 是從屬于 中某向量范數(shù) 的矩陣范數(shù)， 則方程組

3、有唯一不動(dòng)點(diǎn) ，迭代格式 自任意迭代初值 出發(fā)收斂于該不動(dòng)點(diǎn)，且有先驗(yàn)誤差估計(jì)及后驗(yàn)誤差估計(jì)注： q (局部考察q*)越小，迭代收斂越快。11加速迭代法Aitken 加速迭代公式 方程組的加速迭代公式： k = 0, 1, 2, . . 標(biāo)量方程的加速迭代公式：12Newton迭代法f (x) = 0 Newton迭代法： 迭代終止準(zhǔn)則： （預(yù)設(shè)精度）或 （預(yù)設(shè)精度） 非線(xiàn)性方程組的Newton迭代法：F (X) = 013簡(jiǎn)化Newton法 & Newton下山法 簡(jiǎn)化的Newton迭代法 Newton下山法基本思想：減少導(dǎo)數(shù)的計(jì)算 下山因子的取法： 從 =1 開(kāi)始，逐次減半，直到滿(mǎn)足下山條

4、件 基本思想：要求每一步迭代滿(mǎn)足下山條件 具體做法：加下山因子保證全局收斂14簡(jiǎn)化Newton法 & Newton下山法 簡(jiǎn)化的Newton迭代法 Newton下山法基本思想：減少導(dǎo)數(shù)的計(jì)算 下山因子的取法： 從 =1 開(kāi)始，逐次減半，直到滿(mǎn)足下山條件 基本思想：要求每一步迭代滿(mǎn)足下山條件 具體做法：加下山因子保證全局收斂15Jacobi矩陣的計(jì)算 計(jì)算Jacobi矩陣的函數(shù)% 求函數(shù)fun在x處的Jacobi矩陣function J=Jacobi(fun,x)m=length(x); % 確定變量的維度n=length(fun(x); % 確定函數(shù)的維度J=zeros(m,n); % 初始化

5、Jacobi矩陣epsnew=1.d-12; % 小增量I=eye(m,n); % 保存單位向量的矩陣for k=1:n % 依次計(jì)算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)16標(biāo)量方程的Newton迭代算法% fun和dfun以字符串的形式輸入即可function xvect,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,nit=newton(0.5,