課程素材MATLAB在數(shù)學(xué)中的應(yīng)用六章圖

第9章 若有函數(shù)y=f(x)，給定一個(gè)區(qū)間[a,b]n+1個(gè)點(diǎn)x0,x1,…,xna≤x0<x1<…<xn≤b,要求用p(x)f(x)p(xi)=f(xi)=yi,p(x)f(x)的插值多項(xiàng)式。假設(shè)p(x)=a0xn+a1xn-1+…+an-1x+ana0,a1,a2,…,an，由于有p(xi)=f(xi)=yi，于是可以解一個(gè)n+1個(gè)未知數(shù)的線性方程組。…axn+axn- x+a0 1 n-1

xx1x0

xx0x1顯然L1(x0)=f(x0)=y0,L1(x1)=f(x1)=y1,滿足插值條件，所以L1(x)xx1x0

xx0x1

,則稱l0(x)和l1(x)為x0與x1于是有當(dāng)n=23點(diǎn)

(xx1)(xx2)(x0x1)(x0x2

(xx0)(xx2

(x1x0)(x1x2(xx1)(xx0)(x2x1)(x2x0稱為關(guān)于點(diǎn)x0,x1,x2

是，得到的二次插值多項(xiàng)式L2(x)可表示3n次Lagrangen已知n+1個(gè)點(diǎn)(xi,f(xi))(i=0,1,…,n)Ln(x),應(yīng)滿足的條件為Ln(xi)=f(xi)=yi。用插值基函數(shù)方法可得n次Lagrange插值多項(xiàng)式為nLn(x)li(x)f(xi

(xx0)...(xxi1)(xxi1)...(xxn(xix0)...(xixi1)(xixi1)...(xixn)稱為關(guān)于x0,x1,…,xn的n

定理 1n次Lagrange插值多項(xiàng)式Ln(x)是存在且唯一的證明假設(shè)，n次Lagrange插值多項(xiàng)式Ln(x)為a0xn+a1xn-1+…+an-1x+an,已知n+1個(gè)點(diǎn)(xi,f(xi))(i=0,1,…,n)Ln(x),應(yīng)滿足的條件為Ln(xi)=f(xi)=yin+1個(gè)已知條件帶入，則得到一個(gè)有n+1個(gè)方程n+1個(gè)未知數(shù)的線性方程組，由線性代數(shù)的知識可知道，該方程組有唯一解。故Ln(x)是存在且唯一的。9-2設(shè)f(x)∈Cn+1[a,b](表示f(x)在[a,b]n+1階導(dǎo)數(shù))，且節(jié)點(diǎn)f(n1)()x0<x1<…<xn≤b，則插值多項(xiàng)式Ln(x)f(x)Rn(x)=f(x)-中，a<<b,n(x)(xx0)(xx1)...(xxn

(n

n1證明:由插值條件可知Rn(xi)=0(i=0,1,…n)x∈[a,b]Rn(x)=k(x)(x-x0)(x-x1)…(x-xn)=k(x)n1其中，k(x)是依賴x的待定函數(shù)，將x∈[a,b]看做區(qū)間[a,b]上任一固定點(diǎn)，作函數(shù)(t)f(t)Ln(t)k(x)(tx0)(tx1)...(txn顯然(xi0，(i=0,1,…,n),且(x0,它表示(t在[a,b]n+2rolle定理可知(t在[a,b]n+1rolle定理，可知n1(t在[a,b]上至少有一個(gè)零點(diǎn)∈[a,b]，使(n1)()

f(n1)()(n1)!k(x)

f(n1)((n

9-1sin0.32=0.314567,sin0.34=0.333487,sin0.36=0.352274，用線性插值及二次插值計(jì)算sin0.3367的近似值并估計(jì)誤差。解由題意y=f(x)=sin(x)，x0=0.32,y0=0.314567,

xx1x0

xx0x1將x=0.3367|R1(x)|=f’’(x)/2*(x-x0)(x-x1)(xx1)(xx2

(x

)(xx2

(xx1)(xx0=(x0x1)(x0x2

(x1x0)(x1x2

(x2x1)(x2x0x=0.3367)4.nLagrangelagan.m如下(c為插值多項(xiàng)式的系數(shù)function[c,l]=lagan(x,y)fork=1:n+1forif

命令窗口調(diào)用如下（9-1為例9-19-1lagrange從圖可知，求得NewtonNewton其中ai(i=0,…,n)為待定常數(shù)。所有的ai可以根據(jù)插值條件Nn(xi)=f(xi)=yix=x0時(shí)，a0=f(x0)=y0;x=x1時(shí)，a1=f(x1f(x0x1為了求出a2,a3,…,an9-1記f[xm]=f(xm)f

f[xmf[x0x0,xmk-1xmf[x0,…,xk-1,xk]=f[x1,...,xk1,xk]f[x0,x1,...,xk1xk

稱為fx0,x1,….xk-1,xm…其中)=為Rn(xf[x,x0,…,xnn1(xNewton9-2sin0.32=0.314567,sin0.33=0.324043,sin0.34=0.333487,sin0.36=0.352274,Newtonsin0.3367的近似值。解依題意有x0=0.32,x1=0.33,x2=0.34,x3=0.35,x4=0.36,f(x)9-19-1f(x)------Newton將x=0.3367N4(0.3367)=0.330374Newton自定義函數(shù)如下（c表示多項(xiàng)式系數(shù)，d表示均差表functionforfor

命令窗口運(yùn)行結(jié)果如圖9-2所示9-2Newton從圖可知，求得b在工程應(yīng)用中我們經(jīng)常會遇到定積分的計(jì)算ab

f(x)dxF(x)f(x)在[a,b]baf(x)dx=F(b)-F(a)bnf(x)在積分區(qū)間[a,b]上某些點(diǎn)處的函數(shù)值的線性組合naf(x)dxAif(xi)E(fbbxi[a,b]稱為求積節(jié)點(diǎn)，Aixif(x)的具體表baf(x)dxb

b2

(f(a)fE(f)=-1(ba)3fbaf(x)dxb

b6

(f(a)4f

a2

)f(b))

(b

fNewton-Cotes將區(qū)間[a,b]n等分，取等距節(jié)點(diǎn)：a≤x0<x1<…<xn≤b，其中xi=ba(i1)a，i=1,…,n+1n則Newton-Cotes求積可表示bIb

f(x)dx Aif(xi)nn

nn=1,2,3時(shí)，Newton-Cotes求積分別變成了梯形求積，拋物型求積（1/3simpson求積）和3/8simpson求積。Ai分成Ibf(x)dxh(f2f...2f

)

h=(b-a)/n，xi=a+i-1)/h，fi=f(xi)，i=1,2,…,n+1，E(fbah2f若將區(qū)間[a,b]n等分（n為偶數(shù)，對每二個(gè)小區(qū)間用拋物型求積，則得到復(fù)合拋物型求積（復(fù)合1/3simpson求積。Ibf(x)dxh(

4

2

4

...2

4f

)

h

其中，h=(b-a)/nfi=f(a+(i-1)*h)，E(f)=-(b-

9-3分別用梯形求積、1/3simpson求積、n=6的復(fù)合梯形求積、n=61/3simpson求積求定積分1(x22x1)dx0解由高等數(shù)學(xué)知識1(x22x1)dx0用梯形求積1(x22x1)dx=(1-01/3simpson求積1(x22x1)dx0n=6012345671(x22x1)dx=1/12*(f+2f+2f+2f+2f+2f+f01234567n=61/3simpson012345671(x22x1)dx=1/18*(f+4f+2f+4f+2f+4f+f01234567求積算法實(shí)9-3自定義函數(shù)fff.m如下functionz=fff(x)functionI=fhtx(a,b,n)fori=1:n+1則得到結(jié)果則得到結(jié)果I=2.5fff的內(nèi)容改掉即可，若要變換區(qū)間，將函數(shù)調(diào)用中第1、2個(gè)參數(shù)改掉即可。1/3simpsonfunctionI=fhsimpson(a,b,n)fori=1:n+1ifelseifelseifmod(i,2)==0,I=I+4*fff(x(i));elseI=I+2*fff(x(i));則得到結(jié)果則得到結(jié)果I=2.33331/3simpsonf(x)x0處進(jìn)行Taylor(xx (xxx0)f’(x0)+ f''(x

)... f(n)(x

)k導(dǎo)數(shù)定義為fx)k

f(xkh)f(xkh'則得到兩 的數(shù)值微分fk'

fk1h

O(hO(h)f(xk1)f(xk)h

f'

)

f''

的

fxkkkkkkf(xk1)f(xk)

f'(x)

h...

f

)kk'則得到三點(diǎn)的數(shù)值微分fk'

fk1fk1O(h2)kf'k

(fk

8

k

8

k

fk

)O(h4三點(diǎn)：

''

fk12fkfk1O(h2五點(diǎn)：f

''

(

k

16

k

30

16

k

fk

)O(h4x*f(x*)=0，則稱x*f(x)=0f(x)的零點(diǎn)，如果函數(shù)f(x)其中g(shù)(x*)≠0，m為正整數(shù)，當(dāng)m=1時(shí)，稱x*為f(x)=0的單根或f(x)的單零點(diǎn)；當(dāng)m≥2時(shí)，稱x*為f(x)=0的m重根或f(x)的m重零點(diǎn)。如果f(x)為函數(shù)，則稱f(x)=0為方程；如果f(x)為n次多項(xiàng)式，則稱f(x)=0為n次代數(shù)方程。3x-cosx-1=03x-1=cosx,y=3x-1y=cosx的x*∈[1/3,1]。3x-cosx+0.01sinx-1=03x-cosx-1=0在某一區(qū)間上以適當(dāng)?shù)牟介Lh,去函數(shù)值f(xi)(xi=x0+ih,i=0,1,2,…)的符號，當(dāng)f(x)=0。設(shè)函數(shù)f(x)∈C[a,b]二分法的基本思想是：用對分區(qū)間的方法根據(jù)分點(diǎn)處函數(shù)f(x)的符號逐步將有根區(qū)間為方便起見，記a0=a,b0=b。用中點(diǎn)x0=(a0+b0)/2將區(qū)間[a0,b0]分成2個(gè)小區(qū)間[a0,x0]和[x0,b0]。計(jì)算f(x0)。若f(x0)=0，則x0為f(x)=0的根，計(jì)算結(jié)束；否則若f(a0)f(x0)<0，令a1=a0,b1=x0;若f(b0)f(x0)<0，令a1=x0,b1=b0,不管那種情況，都有 bn-an2(bn1an122(bn2an22n(ba

(bafunction[c,err,yc]=bisect(f,a,b,delta)iffork=1:max1ifyc==0elseifyb*yc>0ifb-9-4求f(x)=x3-x-1=0在[1，1.5]內(nèi)的根，要求精確度達(dá)到10-4。可以先定義函數(shù)f如下:functiony=f(x)i為求得的根，j為最后一個(gè)區(qū)間的長度，k為所求點(diǎn)的函數(shù)值（0若有方程f(x)=0在區(qū)間[a,b]1x*f(x)=0改寫成等價(jià)的形式x=(x)，取x0∈[a,b]，利用遞推xk+1=(xk)可以求得序列x1,x2,…,xn。如果迭代方法是收斂的，則xn就是所求的根。k9-5f(x)=x3-x-1=0x0=1.5]附近的根。解若將方程改為x=x3-1,構(gòu)造迭代kxk

x3由x0=1.5,通過迭代可得到x1=2.375,x2=12.39,…,可知，xk顯然不收斂。若將方程改為x=(x+1)1/3,構(gòu)造迭代9-3設(shè)C[a,b],如果對x[a,b]a(x)≤bL(0,1)使|(y)|≤L|x-y|,x,y[a,b]，則在區(qū)間[a,b]上存在惟一點(diǎn)x*,使生成的迭代序列{xk}*x0[a,b]x*，并有誤差估計(jì)|xk*

|1L|x1x0證明先證明x*的存在性，記f(x)=x-(x)f(a)=a-(a)≤0f(b)=b-(b)≥0f(a)=0f(b)=0x*x* |x*x*||(x*)(x*)|L|x*x*||x*x*|與假設(shè)，所以x*x*，即 kx*kkxk

k

L

xx*0因?yàn)?<L<1, lim|xk-x*|=0，即lim0

1=1

|xk1xk|1L|x1x0推論若C1[a,b]（表示在[a,b]上一階導(dǎo)數(shù)連續(xù)，則定理9-3max|(x|L110.3039-6構(gòu)造不同迭代法求x2-3=0的根3解(1)

k13xk3x

(k=0,1,…),(x)=3，'(x)

,(x*1(2)

x1(x2

k

(x)x1(x23),'(x)11x'

31

4xk

1

23x

2 k(x)1(x3),'(x)1(13),'(x*)'

3)01 x29-6m文件如下functiony1=g1(x,n)fori=1:n在命令窗口中調(diào)用的m文件如下functiony2=g2(x,n)fori=1:n 1.7321。的m文件如下functiony3=g3(x,n)fori=1:n NewtonNewton在xkf(x)

f(xk)

(xk)(xxk)

f''()

(xxk

'2kx'2k

f(xk

f'(x

xk+1，于是得到法迭代

x

f(xkk

k f'(xkfunction%ff(x),df表示函數(shù)f'(x)，p0表示迭代的初值，epsilonforf(x)=x3-x-1=01.5附近的根，可以先定義函數(shù)f(x)f’(x)functiony=f(x)functiony=df(x)則得到的結(jié)果z=1.3247f(x)GaussGauss0………function%a為系數(shù)矩陣，b為常量項(xiàng)，n為方程的數(shù)目,相當(dāng)于線性方程AX=bfork=1:nforj=n:-fori=k+1forj=n:-fork=n-1:-1:1fora=[234;352;4330]最后得到結(jié)果x=[-138列主元Gauss列主元Gauss能變?yōu)?，這時(shí)，Gauss消去法無法完成，于是可以考慮列主元Gauss消去法。GaussiiiGauss消去法消元。對每一行都重復(fù)具體方法可以將系數(shù)矩陣a與常量項(xiàng)b組成一個(gè)增廣矩陣,203142541142542031011217121706 11106014 001X3=-列主元Gaussfunctionx=colgauss(a,b,n)fork=1:nfori=k+1nif

forj=n:-1:kforforj=n:-1:k

fork=n-1:-1:1for9-3guass3.無回帶過程的列主元Gauss functionx=c84(a,b,n)fork=1:nfori=k+1nifforj=n:-

forforj=n:-1:k

fori=n:-fork=i-1:- 9-4Gauss4.LULU若將系數(shù)矩陣A分解成單位下三角矩陣L和上三角矩陣Ax=bLUx=b，這時(shí)，可以分解成解方程組Ly=b和Ux=y?？衫萌缦虑蟪鯨和U矩陣

,li1

,i

r lrkuki,r2,3,...,n,ir,rkr

(airlikukr)/urr,r2,3,...,n,ir,rkLy=byk

k lkryr,r

最后利用Ux=yk 1(kukk

nnk

r

LU分解的算法實(shí)現(xiàn)function[x]=LU1(a,b,n)fori=1:nforfori=r

forforr=1:k-

fork=n:-1:1for

Jacobi設(shè)有n…若系數(shù)矩陣非奇異，且aii≠0i=1,2,…,n)

2x=2

(b-

x-

x-…-ax

21

23

2n…nx=1(b-n

x-

x-…-

n1

n2

n-1x(k1)=1

(b1-a12x(k)-…-a1nx(k)2n2n2ax(k1)=2a

(b2-a21x(k)-a23x(k)-…-a2nx(k) …nx(k1)n

1(b-

x(k)

n1

n2

寫成向量形式為其中，B=I-D-1A,F=D-1b,D為矩陣A的對角線元素所組成的對角矩陣，該迭代稱為Jacobi迭function%a為系數(shù)矩陣，b為常量項(xiàng)，n為方程的數(shù)目，m為最大的迭代次數(shù),xforforforif10x1-x2--x1+10x2- 9-5jacobi從圖可知，得到方程組的解為在Jacobi迭代中如果迭代是收斂的，x(k1)應(yīng)該比x(k)更接近方程的解x（i=1,2,…,n 因此，在迭代過程中及時(shí)地用已經(jīng)求出的xk1來代替xki=1,2,…,n-1) 于是可以將Jacobi1x(k1)=1

(b1-a12x(k)-…-a1nx(k)2n2n2ax(k1)=2a

(b2-a21x(k1)-a23x(k)-…-a2nx(k) …nxkn

1(b-

xk1)(-

xk

n1

n2

該迭代稱為Gauss-Seidelixk1)(1i

ii

k1)(

k

ijj

ijnn前面介紹的Jacobi迭成矩陣形式為x(k+1)=Bx(k)+F，Gauss-Seidel迭代也可以寫成矩陣形式x(k+1)=(I-L)-1Ux(k)+(I-L)-1F，其中B=L+U，L為下三角矩陣，U是上三角矩陣。functiona為系數(shù)矩陣，b為常量項(xiàng)，n為方程的數(shù)目，m為最多的迭代次數(shù)，xforfor

forj=i+1 9-6Gauss-Seidel得到方程組的解為i在Gauss-Seidelx(k1)1i

x(k1)

ax(k)),i~(k

n

ijj

ijnjnjj

aijj1

x(k1)jjj

x(k)),i將迭代進(jìn)行加速x(k1)1)x(k)

~(k

ix(k1)(1)x(k)

i

nx(k1)n

ax(k)),ia a

ijj

ijjX(k+1)=Bx(k)+FB=(IL)-11-)I+UF(IL)-1b。若=1，該迭代變成了Gauss-Seidel迭代。functionw為松弛因子，其它與Gauss-Seidelforforfor給定松弛因子w=1.45， 9-7得到方程組的解為dyy'

f(x,y),且

)

（已知

)f(x,

)y(x1)y(x0

，令h=x-x

1 y(x1)=y(x0)+hf(x0,y0)y1f(x1,y1)y(x2)=y(x1)+hf(x1,y1)，一般地，有y(xn+1)=y(xn)+hf(xn,yn)，n=0,1,…,N-1，h=xn+1-xn?？梢杂泍(xn+1yn+1，則有yn+1=yn+hf(xn,yn)Euler剛才討論的就是Euler方%格式 dyfun為函數(shù)f(x,y),xspan為求解%間[x0,xN],y0為初值y(x0),h為步長，x返回節(jié)點(diǎn),yforn=1:length(x)-改進(jìn)的Euleryn+1=yn+h/2[f(xn,ynf(xn+1,yn+1)]，n=0,1,2,…,N-1。%用途：改進(jìn)Euler格式微分方%格式 dyfun為函數(shù)f(x,y),xspan為求解區(qū) [x0,xN],y0為初值y(x0),h為步長，x返回節(jié)點(diǎn),yforn=1:length(x)-k2=feval(dyfun,x(n+1),y(n+1));隱式Euler%用途：隱式Euler格式微分方%格式 dyfun為函數(shù)f(x,y),xspan為求解區(qū) [x0,xN],y0為初值y(x0),h為步長，x返回節(jié)點(diǎn),yforn=1:length(x)-whileabs(y-y1)>ey=y0+k=k+1 變形Euleryn+1=yn+hf(xn+h/2,yn+h/2f(xn,yn))，n=0,1,2,…,N-1。forn=1:length(x)-1Heunyn+1=yn+h/4[f(xn,yn)+3f(xn+2/3h,yn+2/3hf(xn,yn))]，n=0,1,2,…,N-1。forn=1:length(x)-

經(jīng)典四階Runge-Kutta1212

1212

12

forn=1:length(x)-1y’=y-y在[0，1]中的數(shù)值解（取h=0.2Euler方法、改進(jìn)的Euler方法、隱式Euler方法、變形的Euler方法、Heun方法、經(jīng)典四階Runge-Kutta方法等方法求解，并與真解進(jìn)