课程素材MATLAB在数学中的应用六章图

第9章 若有函数y=f(x)，给定一个区间[a,b]n+1个点x0,x1,…,xna≤x0<x1<…<xn≤b,要求用p(x)f(x)p(xi)=f(xi)=yi,p(x)f(x)的插值多项式。假设p(x)=a0xn+a1xn-1+…+an-1x+ana0,a1,a2,…,an，由于有p(xi)=f(xi)=yi，于是可以解一个n+1个未知数的线性方程组。…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于是有当n=23点

(xx1)(xx2)(x0x1)(x0x2

(xx0)(xx2

(x1x0)(x1x2(xx1)(xx0)(x2x1)(x2x0称为关于点x0,x1,x2

是，得到的二次插值多项式L2(x)可表示3n次Lagrangen已知n+1个点(xi,f(xi))(i=0,1,…,n)Ln(x),应满足的条件为Ln(xi)=f(xi)=yi。用插值基函数方法可得n次Lagrange插值多项式为nLn(x)li(x)f(xi

(xx0)...(xxi1)(xxi1)...(xxn(xix0)...(xixi1)(xixi1)...(xixn)称为关于x0,x1,…,xn的n

定理 1n次Lagrange插值多项式Ln(x)是存在且唯一的证明假设，n次Lagrange插值多项式Ln(x)为a0xn+a1xn-1+…+an-1x+an,已知n+1个点(xi,f(xi))(i=0,1,…,n)Ln(x),应满足的条件为Ln(xi)=f(xi)=yin+1个已知条件带入，则得到一个有n+1个方程n+1个未知数的线性方程组，由线性代数的知识可知道，该方程组有唯一解。故Ln(x)是存在且唯一的。9-2设f(x)∈Cn+1[a,b](表示f(x)在[a,b]n+1阶导数)，且节点f(n1)()x0<x1<…<xn≤b，则插值多项式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的待定函数，将x∈[a,b]看做区间[a,b]上任一固定点，作函数(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]上至少有一个零点∈[a,b]，使(n1)()

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

f(n1)((n

9-1sin0.32=0.314567,sin0.34=0.333487,sin0.36=0.352274，用线性插值及二次插值计算sin0.3367的近似值并估计误差。解由题意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为插值多项式的系数function[c,l]=lagan(x,y)fork=1:n+1forif

命令窗口调用如下（9-1为例9-19-1lagrange从图可知，求得NewtonNewton其中ai(i=0,…,n)为待定常数。所有的ai可以根据插值条件Nn(xi)=f(xi)=yix=x0时，a0=f(x0)=y0;x=x1时，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自定义函数如下（c表示多项式系数，d表示均差表functionforfor

命令窗口运行结果如图9-2所示9-2Newton从图可知，求得b在工程应用中我们经常会遇到定积分的计算ab

f(x)dxF(x)f(x)在[a,b]baf(x)dx=F(b)-F(a)bnf(x)在积分区间[a,b]上某些点处的函数值的线性组合naf(x)dxAif(xi)E(fbbxi[a,b]称为求积节点，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将区间[a,b]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时，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若将区间[a,b]n等分（n为偶数，对每二个小区间用抛物型求积，则得到复合抛物型求积（复合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的复合梯形求积、n=61/3simpson求积求定积分1(x22x1)dx0解由高等数学知识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求积算法实9-3自定义函数fff.m如下functionz=fff(x)functionI=fhtx(a,b,n)fori=1:n+1则得到结果则得到结果I=2.5fff的内容改掉即可，若要变换区间，将函数调用中第1、2个参数改掉即可。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));则得到结果则得到结果I=2.33331/3simpsonf(x)x0处进行Taylor(xx (xxx0)f’(x0)+ f''(x

)... f(n)(x

)k导数定义为fx)k

f(xkh)f(xkh'则得到两 的数值微分fk'

fk1h

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

f'

)

f''

的

fxkkkkkkf(xk1)f(xk)

f'(x)

h...

f

)kk'则得到三点的数值微分fk'

fk1fk1O(h2)kf'k

(fk

8

k

8

k

fk

)O(h4三点：

''

fk12fkfk1O(h2五点：f

''

(

k

16

k

30

16

k

fk

)O(h4x*f(x*)=0，则称x*f(x)=0f(x)的零点，如果函数f(x)其中g(x*)≠0，m为正整数，当m=1时，称x*为f(x)=0的单根或f(x)的单零点；当m≥2时，称x*为f(x)=0的m重根或f(x)的m重零点。如果f(x)为函数，则称f(x)=0为方程；如果f(x)为n次多项式，则称f(x)=0为n次代数方程。3x-cosx-1=03x-1=cosx,y=3x-1y=cosx的x*∈[1/3,1]。3x-cosx+0.01sinx-1=03x-cosx-1=0在某一区间上以适当的步长h,去函数值f(xi)(xi=x0+ih,i=0,1,2,…)的符号，当f(x)=0。设函数f(x)∈C[a,b]二分法的基本思想是：用对分区间的方法根据分点处函数f(x)的符号逐步将有根区间为方便起见，记a0=a,b0=b。用中点x0=(a0+b0)/2将区间[a0,b0]分成2个小区间[a0,x0]和[x0,b0]。计算f(x0)。若f(x0)=0，则x0为f(x)=0的根，计算结束；否则若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]内的根，要求精确度达到10-4。可以先定义函数f如下:functiony=f(x)i为求得的根，j为最后一个区间的长度，k为所求点的函数值（0若有方程f(x)=0在区间[a,b]1x*f(x)=0改写成等价的形式x=(x)，取x0∈[a,b]，利用递推xk+1=(xk)可以求得序列x1,x2,…,xn。如果迭代方法是收敛的，则xn就是所求的根。k9-5f(x)=x3-x-1=0x0=1.5]附近的根。解若将方程改为x=x3-1,构造迭代kxk

x3由x0=1.5,通过迭代可得到x1=2.375,x2=12.39,…,可知，xk显然不收敛。若将方程改为x=(x+1)1/3,构造迭代9-3设C[a,b],如果对x[a,b]a(x)≤bL(0,1)使|(y)|≤L|x-y|,x,y[a,b]，则在区间[a,b]上存在惟一点x*,使生成的迭代序列{xk}*x0[a,b]x*，并有误差估计|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*|与假设，所以x*x*，即 kx*kkxk

k

L

xx*0因为0<L<1, lim|xk-x*|=0，即lim0

1=1

|xk1xk|1L|x1x0推论若C1[a,b]（表示在[a,b]上一阶导数连续，则定理9-3max|(x|L110.3039-6构造不同迭代法求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在命令窗口中调用的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表示函数f'(x)，p0表示迭代的初值，epsilonforf(x)=x3-x-1=01.5附近的根，可以先定义函数f(x)f’(x)functiony=f(x)functiony=df(x)则得到的结果z=1.3247f(x)GaussGauss0………function%a为系数矩阵，b为常量项，n为方程的数目,相当于线性方程AX=bfork=1:nforj=n:-fori=k+1forj=n:-fork=n-1:-1:1fora=[234;352;4330]最后得到结果x=[-138列主元Gauss列主元Gauss能变为0，这时，Gauss消去法无法完成，于是可以考虑列主元Gauss消去法。GaussiiiGauss消去法消元。对每一行都重复具体方法可以将系数矩阵a与常量项b组成一个增广矩阵,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若将系数矩阵A分解成单位下三角矩阵L和上三角矩阵Ax=bLUx=b，这时，可以分解成解方程组Ly=b和Ux=y。可利用如下求出L和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分解的算法实现function[x]=LU1(a,b,n)fori=1:nforfori=r

forforr=1:k-

fork=n:-1:1for

Jacobi设有n…若系数矩阵非奇异，且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为系数矩阵，b为常量项，n为方程的数目，m为最大的迭代次数,xforforforif10x1-x2--x1+10x2- 9-5jacobi从图可知，得到方程组的解为在Jacobi迭代中如果迭代是收敛的，x(k1)应该比x(k)更接近方程的解x（i=1,2,…,n 因此，在迭代过程中及时地用已经求出的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为系数矩阵，b为常量项，n为方程的数目，m为最多的迭代次数，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将迭代进行加速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。可以记y(xn+1yn+1，则有yn+1=yn+hf(xn,yn)Euler刚才讨论的就是Euler方%格式 dyfun为函数f(x,y),xspan为求解%间[x0,xN],y0为初值y(x0),h为步长，x返回节点,yforn=1:length(x)-改进的Euleryn+1=yn+h/2[f(xn,ynf(xn+1,yn+1)]，n=0,1,2,…,N-1。%用途：改进Euler格式微分方%格式 dyfun为函数f(x,y),xspan为求解区 [x0,xN],y0为初值y(x0),h为步长，x返回节点,yforn=1:length(x)-k2=feval(dyfun,x(n+1),y(n+1));隐式Euler%用途：隐式Euler格式微分方%格式 dyfun为函数f(x,y),xspan为求解区 [x0,xN],y0为初值y(x0),h为步长，x返回节点,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)-

经典四阶Runge-Kutta1212

1212

12

forn=1:length(x)-1y’=y-y在[0，1]中的数值解（取h=0.2Euler方法、改进的Euler方法、隐式Euler方法、变形的Euler方法、Heun方法、经典四阶Runge-Kutta方法等方法求解，并与真解进