feval

1

非线性方程数值方法

1、固定点迭代：求解方程)(xgx的近似值，初始值为

0

x，迭代式为)(

1nn

xgx



function[k,x1,err,x]=fixpt(g,x0,tol,maxl)

%Input-gistheiterationfunctioninputasastring'g'

%-x0istheinitialguessforthefixedpoint

%-tolisthetolerance

%-maxlisthemaximumnumberofiterations

%Output-kisthenumberofiterationthatwerecarriedout

%-x1istheapproximationtothefixedpoint

%-erristheerrorintheapproximation

%-xcontainsthequence{xn}

x(1)=x0;

fork=2:maxl

x(k)=feval(g,x(k-1));

err=abs(x(k)-x(k-1));

relerr=err/(abs(x(k))+eps);

x1=x(k);

if(err

end

ifk==maxl

disp('maximumnumberofiterationxceeded')

end

x=x';

例如g=inline('exp(-x)');

[k,p,err,P]=fixpt(g,0.5,0.01,20)

2、二分法：求解方程0)(xf在区间],[ba内的一个根。前提条件是)(xf是连续的，且)(af

与)(bf的符号相反。

function[c,err,yc]=bict(f,a,b,delta)

%Input-fisthefunctioninputasastring'f'

%-aandbaretheleftandrightendpoint

%-deltaisthetolerance

%Output-cisthezero

%-yc=f(c)

%-erristheerrorestimateforc

ya=feval(f,a),yb=feval(f,b)

ifya*yb>0,break,end

maxl=1+round((log(b-a)-log(delta))/log(2))

fork=1:maxl

c=(a+b)/2;yc=feval(f,c);

ifyc==0

a=c;b=c;

elifyb*yc>0

b=c;yb=yc;

el

a=c;ya=yc;

end

ifb-a

end

c=(a+b)/2;err=abs(b-a);yc=feval(f,c);

例：f=inline('x*sin(x)-1')

2

[c,err,yc]=bict(f,0,2,0.01)

3、牛顿－拉夫申迭代：使用初始近似值

0

x，利用迭代式,2,1,

)('

)(

1

1

1









k

xf

xf

xx

k

k

kk

，计算

函数0)(xf的根的近似值。

function[x0,err,k,y]=newton(f,df,x0,delta,epsilon,maxl)

%Input-fistheobjectfunctioninputasastring'f'

%-dfisthederivativeoffinputasastring'df'

%-x0istheinitialapproximationtoazerooff

%-deltaisthetoleranceforx0

%-epsilonisthetoleranceforthefunctionvaluesy

%-maxlisthemaximumnumberofiterations

%Output-xistheNewton-Raphsonapproximationtothezero

%-erristheerrorestimateforx0

%-kisthenumberofiterations

%-yisthefunctionvaluef(x0)

fork=1:maxl

x1=x0-feval(f,x0)/feval(df,x0);

err=abs(x1-x0);

relerr=2*err/(abs(x1)+delta);

x0=x1;

y=feval(f,x0);

if(err

end

例：f=inline('x^2-5'),df=inline('2*x')

[x,err,k,y]=newton(f,df,2,0.001,0.001,20)

4、割线法：使用初始近似值

0

x和

1

x，利用迭代式,2,1,

)()(

))((

1

1

1















k

xfxf

xxxf

xx

kk

kkk

kk

，计

算函数0)(xf的根的近似值。

function[x1,err,k,y]=cant(f,x0,x1,delta,epsilon,maxl)

%Input-fistheobjectfunctioninputasastring'f'

%-x0andx1istheinitialapproximationstoazero

%-deltaisthetoleranceforx1

%-epsilonisthetoleranceforthefunctionvaluesy

%-maxlisthemaximumnumberofiterations

%Output-x1isthecantmethodapproximationtothezero

%-erristheerrorestimateforx1

%-kisthenumberofiterations

%-yisthefunctionvaluef(x1)

fork=1:maxl

x2=x1-feval(f,x1)*(x1-x0)/(feval(f,x1)-feval(f,x0));

err=abs(x2-x1);

relerr=2*err/(abs(x1)+delta);

x0=x1;x1=x2;

y=feval(f,x1);

if(err

end

例：f=inline('x^3-3*x+2')

[x1,err,k,y]=cant(f,-2.6,-2.4,0.01,0.001,20)

3

5、Steffenn加速法：给定初始近似值

0

p，快速寻找固定点方程)(xgx的解，假设)(xg和

)('xg是连续的，而且1)('xg，通常的固定点迭代缓慢（线性）收敛到p。

function[P,Q]=steff(f,df,p0,delta,epsilon,maxl)

%Input-fistheobjectfunctioninputasastring'f'

%-dfisthederivativeoffinputasastring'df'

%-p0istheinitialapproximationtoazerooff

%-deltaisthetoleranceforp0

%-epsilonisthetoleranceforthefunctionvaluesy

%-maxlisthemaximumnumberofiterations

%Output-PistheSteffennapproximationtothezero

%-QisthematrixcontainingtheSteffennquence

%InitializethematrixR

R=zeros(maxl,3);

R(1,1)=p0;

fork=1:maxl

forj=2:3

%DenominatorinNewton-Raphsonmethodiscalculated

nrodenom=feval(df,R(k,j-1));

%CalculateNewton-Raphsonapproximations

ifnrodenom==0

'divisionbyzeroinNewton-Raphsonmethod'

break

el

R(k,j)=R(k,j-1)-feval(f,R(k,j-1))/nrodenom;

end

%DenominatorinAitken'sAccelerationprocesscalculated

aadenom=R(k,3)-2*R(k,2)+R(k,1);

%CalculateAitken'sAccelerationapproximations

ifaadenom==0

'divisionbyzeroinAitken''sAcceleration'

break

el

R(k+1,1)=R(k,1)-(R(k,2)-R(k,1))^2/aadenom;

end

end

%Endprogramifdivisionbyzerooccurred

if(nrodenom==0)|(aadenom==0)

break

end

%Stoppingcriteriaareevluated

err=abs(R(k,1)-R(k+1,1));

relerr=err/(abs(R(k+1,1))+delta);

y=feval(f,R(k+1,1));

if(err

%PandthematrixQaredetermined

P=R(k+1,1);Q=R(1:k+1,:);

break

end

end

end

例：f=inline('x^2-5'),df=inline('2*x')

[P,Q]=steff(f,df,2,0.001,0.001,20)

4

6、Muller法：给定三个初始值

210

,,ppp，求方程0)(xf的根。

Muller法是一种迭代方法，构造一条抛物线经过初始三点，然后利用二次函数求下一个近似值的二

次根。可以证明当接近一个单根时，Muller法比割线法收敛得快，基本上与牛顿法一样快。此方法

可用来求函数的实数零点和复数零点，且可用于为复杂的算式编程序。

设

2

p是根的最佳近似值。改变变量

2

pxt。使用差分

200

pph和

211

pph。设

包含t的二次多项式为：cbtaty2

。根据三点0,,

10

hht可得到三个包含cba,,的三个方

程，解得

2

fc，并利用定义cfe

00

和cfe

11

，可得

2

10

2

01

2

10

2

01

2

10

2

01

0110,

hhhh

hehe

b

hhhh

hehe

a











（1）

二次式的根为

acbb

c

z

4

2

2



（2）

为确保方法的稳定性，需选择式是绝对值最小的根。如果0b，使用带正号的根；如果0b，

使用带负号的根。则最新的根为zpp

23

。

为更新迭代，需要从},,{

210

ppp中选择最靠近

3

p的两点为新的

10

,pp（即放弃离

3

p最远的

一点）。

function[p,y,err]=muller(f,p0,p1,p2,delta,epsilon,maxl)

%Input-fistheobjectfunctioninputasastring'f'

%-p0,p1,andp2aretheinitialapproximations

%-deltaisthetoleranceforp0,p1,andp2

%-epsilonisthetoleranceforthdfunctionvaluesy

%-maxlisthemaximumnumberofiterations

%Output-pistheMullerapproximationtothezerooff

%-yisthefunctionvaluey=f(p)

%-erristheerrorintheapproximationofp

%InitializethematrixesPandY

P=[p0p1p2];

Y=feval(f,P);

%Calculateaandbinformula(1)

fork=1:maxl

h0=P(1)-P(3);h1=P(2)-P(3);e0=Y(1)-Y(3);e1=Y(2)-Y(3);c=Y(3);

denom=h1*h0^2-h0*h1^2;

a=(e0*h1-e1*h0)/denom;

b=(e1*h0^2-e0*h1^2)/denom;

%Suppressanycomplexroots

ifb^2-4*a*c>0

disc=sqrt(b^2-4*a*c);

el

disc=0;

end

%Findthesmallestroot(2)

ifb<0

disc=-disc;

end

z=-2*c/(b+disc);

p=P(3)+z;

%sorttheentriesofPtofindthetwoclosttop

ifabs(p-P(2))

Q=[P(2)P(1)P(3)];P=Q;Y=feval(f,P);

5

end

%ReplacetheentryofPthatwasfarthestfromPwithp

P(3)=p;Y(3)=feval(f,P(3));

y=Y(3);

%Determinestoppingcriteria

err=abs(z);

relerr=err/(abs(p)+delta);

if(err

break

end

end

例：f=inline('x.^3-3.*x+2');

[p,y,err]=muller(f,1.4,1.3,1.2,0.01,0.001,20)

线性方程组数值方法

1、回代：用回代法求解上三角线性方程组bAx，必须满足系数矩阵的对角元素非零的条

件。

functionx=backsub(A,b)

%Input-Aisannxnupper-triangularnonsingularmatrix

%-bisannx1matrix

%Output-xisthesolutiontothelinearsystemAx=b

%Findthedimensionofbandinitializex

n=length(b);x=zeros(n,1);

x(n)=b(n)/A(n,n);

fork=n-1:-1:1

x(k)=(b(k)-A(k,k+1:n)*x(k+1:n))/A(k,k);

end

例：A=[1214;0-42-5;00-5-7.5;000-9];b=[13;2;-35;-18];

x=backsub(A,b)

2、高斯消去法(不选主元素)

functionx=gauss(A,b)

%Input-Aisannxnnonsingularmatrix

%-bisannx1matrix

%Output-xisannx1matrixcontainingthesolutiontoAx=b

%Initializex

[n,n]=size(A);x=zeros(n,1);

%Formtheaugmentedmatrix:Aug=[A|b]

Aug=[Ab];

fork=1:n-1

%Eliminationprocessforcolumnk

fori=k+1:n

m=Aug(i,k)/Aug(k,k);

Aug(i,k:n+1)=Aug(i,k:n+1)-m*Aug(k,k:n+1);

end

end

%OutputEliminationmatrix

Aug

%BackSubstitution

A=Aug(1:n,1:n);b=Aug(1:n,n+1);

6

x(n)=b(n)/A(n,n);

fork=n-1:-1:1

x(k)=(b(k)-A(k,k+1:n)*x(k+1:n))/A(k,k);

end

例：A=[4-92;2-46;-11-3];b=[5;3;-4];x=gauss(A,b)

3、列主元消去法

functionx=gauss_column(A,b)

%Input-Aisannxnnonsingularmatrix

%-bisannx1matrix

%Output-xisannx1matrixcontainingthesolutiontoAx=b

%InitializexandthetemporarystoragematrixC

[n,n]=size(A);x=zeros(n,1);C=zeros(1,n+1);

%Formtheaugmentedmatrix:Aug=[A|b]

Aug=[Ab];

fork=1:n-1

%Partialpivotingforcolumnk

[y,j]=max(abs(Aug(k:n,k)));

%Interchangerowkandj

C=Aug(k,:);Aug(k,:)=Aug(j+k-1,:);Aug(j+k-1,:)=C;

ifAug(k,k)==0

'Awassingular,Nouniquesolution'

break

end

%Eliminationprocessforcolumnk

fori=k+1:n

m=Aug(i,k)/Aug(k,k);

Aug(i,k:n+1)=Aug(i,k:n+1)-m*Aug(k,k:n+1);

end

end

%BackSubstitutionon[U|y]usingbacksub

x=backsub(Aug(1:n,1:n),Aug(1:n,n+1));

例：A=[1214;2043;4221;-3132];b=[13;28;20;6];

x=gauss_column(A,b)

4、LU分解法(不选主元素)

functionx=lufact(A,b)

%Input-Aisannxnmatrix

%-bisannx1matrix

%Output-xisannx1matrixcontainingthesolutiontoAx=b

%Initilizex,y,

[n,n]=size(A);x=zeros(n,1);y=zeros(n,1);

fork=1:n-1

%CalculatemultiplierandplaceinsubdiagonalportionofA

fori=k+1:n

mult=A(i,k)/A(k,k);

A(i,k)=mult;

A(i,k+1:n)=A(i,k+1:n)-mult*A(k,k+1:n);

end

end

%OutputLUfactorizationcompactmatrix

A

%solvefory

7

y(1)=b(1);

fori=2:n

y(i)=b(i)-A(i,1:i-1)*y(1:i-1);

end

%Solveforx

x(n)=y(n)/A(n,n);

fori=n-1:-1:1

x(i)=(y(i)-A(i,i+1:n)*x(i+1:n))/A(i,i);

end

例：A=[4-92;2-46;-11-3];b=[5;3;-4];x=lufact(A,b)

5、列主元LU分解法（PA＝LU）A是非奇异矩阵

functionx=lu_column(A,b)

%Input-Aisannxnmatrix

%-bisannx1matrix

%Output-xisannx1matrixcontainingthesolutiontoAx=b

%Initilizex,y,thetemporarystoragematrixC,andtherow

%permutationinformationmatrixR

[n,n]=size(A);

x=zeros(n,1);y=zeros(n,1);C=zeros(1,n);R=1:n;

fork=1:n-1

%Findthepivotrowforcolumnk

[maxl,j]=max(abs(A(k:n,k)));

%Interchangerowkandj

C=A(k,:);A(k,:)=A(j+k-1,:);A(j+k-1,:)=C;

d=R(k);R(k)=R(j+k-1);R(j+k-1)=d;

ifA(k,k)==0

'Aissingular,Nouniquesolution'

break

end

%CalculatemultiplierandplaceinsubdiagonalportionofA

fori=k+1:n

mult=A(i,k)/A(k,k);

A(i,k)=mult;

A(i,k+1:n)=A(i,k+1:n)-mult*A(k,k+1:n);

end

end

%OutputLUfactorizationcompactmatrix

A

%solvefory

y(1)=b(R(1));

fori=2:n

y(i)=b(R(i))-A(i,1:i-1)*y(1:i-1);

end

%Solveforx

x(n)=y(n)/A(n,n);

fori=n-1:-1:1

x(i)=(y(i)-A(i,i+1:n)*x(i+1:n))/A(i,i);

end

例：A=[126;48-1;-23-5];b=[1;2;3];

x=lu_column(A,b)

8

6、Jacobi迭代:程序可用的充分条件是A具有严格对角优势。

functionx=jacobi(A,b,x0,delta,maxl)

%Input-Aisannxnnonsingularmatrix

%-bisannx1matrix

%-x0isannx1matrix;theinitialguess

%-deltaisthetoleranceforx0

%-maxlisthemaximumnumberofiterations

%Output-xisannx1matrix;thejacobiapproximationtothe

%solutionofAx=b

n=length(b);

fork=1:maxl

fori=1:n

x(i)=(b(i)-A(i,[1:i-1,i+1:n])*x0([1:i-1,i+1:n]))/A(i,i);

end

err=abs(norm(x'-x0));

relerr=err/(norm(x)+eps);

x0=x';

if(err

break

end

end

x=x';

例：A=[4-11;4-81;-215];b=[7;-21;15];

x=jacobi(A,b,[1;2;2],0.001,50)

7、Gauss-Seidel迭代:程序可用的充分条件是A具有严格对角优势。

functionx=gauss_idel(A,b,x0,delta,maxl)

%Input-Aisannxnnonsingularmatrix

%-bisannx1matrix

%-x0isannx1matrix;theinitialguess

%-deltaisthetoleranceforx0

%-maxlisthemaximumnumberofiterations

%Output-xisannx1matrix;thegauss_idelapproximationtothe

%solutionofAx=b

n=length(b);

fork=1:maxl

fori=1:n

ifi==1

x(1)=(b(1)-A(1,2:n)*x0(2:n))/A(1,1);

elifi==n

x(n)=(b(n)-A(n,1:n-1)*(x(1:n-1))')/A(n,n);

el

%xcontainsthekthapproximationsandx0the(k-1)st

x(i)=(b(i)-A(i,1:i-1)*x(1:i-1)-A(i,i+1:n)*x0(i+1:n))/A(i,i);

end

end

err=abs(norm(x'-x0));

relerr=err/(norm(x)+eps);

x0=x';

if(err

break

end

end

x=x';

9

例：A=[4-11;4-81;-215];b=[7;-21;15];

x=gauss_idel(A,b,[1;2;2],0.001,50)

非线性方程组数值方法

1、非线性Seidel迭代：求解非线性固定点方程组)(xGx。

function[x,iter]=idel(G,x,delta,maxl)

%Input-GisthenonlinearsystemsavedintheM-fileG.m

%-xistheinitialguessatthesolution

%-deltaistheerrorbound

%-maxlisthenumberofiterations

%Output-xistheidelapproximationtothesolution

%-iteristhenumberofiterationsrequired

n=length(x);

fork=1:maxl

x0=x;

%x0isthekthapproximationtothesolution

forj=1:n

A=feval('G',x0);

%Updatethetermsofx0astheyarecalculated

x0(j)=A(j);

end

err=abs(norm(x0-x));

relerr=err/(norm(x0)+eps);

x=x0;iter=k;

if(err

break

end

end

2、牛顿－拉夫申法：求解非线性方程组0)(xF，,2,1,0,

)('

)(

1





k

xF

xF

xx

k

k

kk

。

function[x,iter,err]=newtondim(F,JF,x,delta,epsilon,maxl)

%Input-FisthenonlinearsystemsavedintheM-fileF.m

%-JFistheJacobianofFsavedastheM-fileJF.m

%-xistheinitialapproximationtothesolution

%-deltaisthetoleranceforx

%-epsilonisthetoleranceforF(x)

%-maxlisthenumberofiterations

%Output-xistheidelapproximationtothesolution

%-iteristhenumberofiterationsrequired

%-erristheerrorestimateforx

y=feval(F,x)

fork=1:maxl

J=feval(JF,x);

Q=x-(Jy')';

z=feval(F,Q);

err=norm(Q-x);relerr=err/(norm(Q)+eps);

x=Q;y=z;iter=k;

if(err

break

end

10

end

例：functionz=F(X)

x=X(1);y=X(2);

z=zeros(1,2);

z(1)=x^2-2*x-y+0.5;

z(2)=x^2+4*y^2-4;

functionw=JF(X)

x=X(1);y=X(2);

w=[2*x-2,-1;2*x,8*y];

[x,iter,err]=newtondim('F','JF',[20.25],0.001,0.001,20)

插值方法

1、拉格朗日方法：

function[C,L]=lagrange(x,y)

%Input-xisavectorthatcontainsalistofabscissas

%-yisavectorthatcontainsalistofordinates

%Output-CisamatrixthatcontainsthecoefficientsoftheLagrange

%interpolatorypolynomial

%-LisamatrixthatcontainstheLagrangecoefficientpolynomials

w=length(x);n=w-1;L=zeros(w,w);

%FormtheLagrangecoefficientpolynomials

fori=1:n+1

v=1;

forj=1:n+1

ifi~=j

v=conv(v,poly(x(j)))/(x(i)-x(j));

end

end

L(i,:)=v;

end

%DeterminethecoefficientsoftheLagrangeinterpolatingpolynomial

C=y*L;

例：x=[pi/6pi/4pi/3];y=sin(x);[C,L]=lagrange(x,y)

在z=7*pi/24处的插值为：z=7*pi/24;C*[z^2,z,1]'

2、牛顿插值方法：构造和计算经过点nkyx

kk

,,1,0),,(的次数小于等于n的牛顿插值多项

式：)())(())(()()(

1nnn

xxxxxxdxxxxdxxddxN，

其中

jkk

jkjk

kjkkxx

dd

dyd









1,11,

0

,.

function[C,D]=newtonpoly(x,y)

%Input-xisavectorthatcontainsalistofabscissas

%-yisavectorthatcontainsalistofordinates

%Output-CisamatrixthatcontainsthecoefficientsoftheNewton

%interpolatorypolynomial

%-Disthedivided-differencetable

n=length(x);D=zeros(n,n);D(:,1)=y';

%Formthedivided-differncetable

forj=2:n

11

fork=j:n

D(k,j)=(D(k,j-1)-D(k-1,j-1))/(x(k)-x(k-j+1));

end

end

%DeterminethecoefficientsoftheNewtoninterpolatingpolynomial

C=D(n,n);

fork=(n-1):-1:1

C=conv(C,poly(x(k)));

m=length(C);

C(m)=C(m)+D(k,k);

end

例：x=[pi/6pi/4pi/3];y=sin(x);[C,D]=newtonpoly(x,y)

在z=7*pi/24处的插值为：z=7*pi/24;C*[z^2,z,1]'

3、切比雪夫逼近：构造和计算]1,1[内的N次切比雪夫逼近多项式，其中：







N

j

jj

xTcxP

0

)()(，基于节点：)

22

)12(

cos(







N

k

x

k



。

function[C,X,Y]=chebyshev(fun,n,a,b)

%Input-funisthestringfunctiontobeapproximated

%-nisthedegreeofthechebyshevinterpolatingpolynomial

%-aistheleftendpoint

%-bistherightendpoint

%-Cisthecoefficientlistforthepolynomial

%-Xcontainstheabscissas

%-Ycontainstheordinates

ifnargin==2,a=-1;b=1;end

d=pi/(2*n+2);C=zeros(1,n+1);

fork=1:n+1

X(k)=cos((2*k-1)*d);

end

X=(b-a)*X/2+(a+b)/2;

x=X;

Y=eval('fun');

fork=1:n+1

z=(2*k-1)*d;

forj=1:n+1

C(j)=C(j)+Y(k)*cos(j-1)*z;

end

end

C=2*C/(n+1);

C(1)=C(1)/2;

例：fun=inline('sin(x)');[C,X,Y]=chebyshev(fun,5,0,pi/4)

4、压紧三次样条曲线（带第一类固定边界条件）：根据1n个点

n

kkk

yx

0

),(



，构造并计算三次

样条插值。

12

5、三角多项式逼近（离散傅里叶级数）：设nm12，根据n个等距的横坐标值

nknkx

k

,,2,1,/2，构造m阶三角多项，其形式为：







m

j

jj

jxbjxa

a

xP

1

0))sin()cos((

2

)(，

多项式的系数计算公式为：

mjjxxf

n

a

n

k

kkj

,,1,0,)cos()(

2

1





，mjjxxf

n

b

n

k

kkj

,,1,)sin()(

2

1





。

function[A,B]=tpcoeff(X,Y,m)

%Input-Xisavectorofequallyspacedabscessasin[-pi,pi]

%-Yisavectorofordinates

%-misthedegreeofthetrigonometricpolynomial

%Output-Aisavectorcontainingthecoefficientsofcos(jx)

%-Bisavectorcontainingthecoefficientsofsin(jx)

n=length(X)-1;

maxl=fix((n-1)/2);

ifm>maxl

m=maxl;

end

A=zeros(1,m+1);

B=zeros(1,m+1);

Yends=(Y(1)+Y(n+1))/2;

Y(1)=Yends;

Y(n+1)=Yends;

A(1)=sum(Y);

forj=1:m

A(j+1)=cos(j*X)*Y';

B(j+1)=sin(j*X)*Y';

end

A=2*A/n;

B=2*B/n;

A(1)=A(1)/2;

【例】根据12个等距横坐标点2/)(,12,,2,1,6/xxfkkx

k

，求解5阶三角

多项式逼近。首先利用函数tpcoeff构造m=5阶三角多项式的系数

jj

ba,的矩阵BA,。

X=-pi+[1:12]*pi/6;Y=X/2;n=length(Y);Y(n)=0;

[A,B]=tpcoeff(X,Y,5)

下面程序计算了m=5阶三角多项式)(xP在点x处的值。

functionz=tp(A,B,x,m)

z=A(1);

forj=1:m

z=z+A(j+1)*cos(j*x)+B(j+1)*sin(j*x);

end

在MATLAB命令窗口输入下列命令，画出三角多项式拟合曲线和12个离散点。

x=-pi:.01:pi;y=tp(A,B,x,5);plot(x,y,X,Y,'o')

曲线拟合

1、最小二乘拟合直线：根据n个数据点),(,),,(

11nn

yxyx构造最小二乘拟合直线

baxy。

function[a,b]=lsline(X,Y)

%Input-Xisthe1xnabscissasvector

%-Yisthe1xnordinatevector

13

%Output-aisthecoefficientofxinax+b

%-bistheconstantcoefficientinax+b

xmean=mean(X);

ymean=mean(Y);

sumx2=(X-xmean)*(X-xmean)';

sumxy=(Y-ymean)*(X-xmean)';

a=sumxy/sumx2;

b=ymean-a*xmean;

例：X=[-10123456];Y=[10975430-1];[a,b]=lsline(X,Y)

2、最小二乘多项式拟合：根据n个数据点),(,),,(

11nn

yxyx构造m阶最小二乘多项式

m

m

m

mm

xcxcxcxccxP

1

12

321

)(





functionC=lspoly(X,Y,m)

%Input-Xisthe1xnabscissasvector

%-Yisthe1xnordinatevector

%-misthedegreeoftheleast-squarespolynomial

%Output-Cisthecoefficientlistforthepolynomial

n=length(X);

B=zeros(1:m+1);F=zeros(n:m+1);

%FillthecolumnsofFwiththepowersofX

fork=1:m+1

F(:,k)=X'.^(k-1);

end

%solvethelinearsystem

A=F'*F;

B=F'*Y';

C=AB;

C=flipud(C);

例：X=[-3024];Y=[3113];C=lspoly(X,Y,2)

数值微分

1、使用极限的微分求解：计算)('xf的近似值，生成序列：

nk

h

hxfhxf

Dxf

k

kk

k

,,0,

)10(2

)10()10(

)('









当

11



nnnn

DDDD或

1



nn

DD＜允许误差时停止计算。后一个不等式用来求最侍

近似值

n

Dxf)('。

function[L,n]=difflim(f,x,toler)

%Input-fisthefunctioninputasastring'f'

%-xisthedifferentiationpoint

%-toleristhetolerancefortheerror

%Output-L=[H'D'E']:

%Histhevectorofstepsizes

%Disthevectorofapproximatederivatives

%Eisthevectoroferrorbounds

%-nisthecoordinateofthe"bestapproximation"

maxl=15;

h=1;

14

H(1)=h;

D(1)=(feval(f,x+h)-feval(f,x-h))/(2*h);

E(1)=0;

R(1)=0;

forn=1:2

h=h/10;

H(n+1)=h;

D(n+1)=(feval(f,x+h)-feval(f,x-h))/(2*h);

E(n+1)=abs(D(n+1)-D(n));

R(n+1)=2*E(n+1)*(abs(D(n+1))+abs(D(n))+eps);

end

n=2;

while((E(n)>E(n+1))&(R(n)>toler))&n

h=h/10;

H(n+2)=h;

D(n+2)=(feval(f,x+h)-feval(f,x-h))/(2*h);

E(n+2)=abs(D(n+2)-D(n+1));

R(n+2)=2*E(n+2)*(abs(D(n+2))+abs(D(n+1))+eps);

n=n+1;

end

n=length(D)-1;

L=[H'D'E'];

【例】设

xexf)(，且1x。使用步长10,,2,1,10khk

k

计算差商

k

D。精度为小

数点后9位。

f=inline('exp(x)');x=1;toler=0.5*10^(-9);

[L,n]=difflim(f,x,toler),vpa(L,10)

2、利用Richardson外推法的微分求解：求解)('xf的数值解，构造包含近似值jkkjD),,(

的表，并将),()('nnDxf作为最终答案。近似值),(kjD存放在下三角矩阵中。第一列是：

h

hxfhxf

jD

j

jj

12

)2()2(

)0,(





,

行j的元素为：jk

kjDkjD

kjDkjD

k







1,

14

)1,1()1,(

)1,(),(。

function[D,err,relerr,n]=diffext(f,x,delta,toler)

%Input-fisthefunctioninputasastring'f'

%-xisthedifferentiationpoint

%-deltaisthetolerancefortheerror

%-toleristhetolerancefortherelativeerror

%Output-Disthematrixofapproximatederivatives

%-erristheerrorbound

%-relerristherelativeerrorbound

%-nisthecoordinateofthe"bestapproximation"

err=1;

relerr=1;

h=1;

j=1;

D(1,1)=(feval(f,x+h)-feval(f,x-h))/(2*h);

whilerelerr>toler&err>delta&j<12

h=h/2;

D(j+1,1)=(feval(f,x+h)-feval(f,x-h))/(2*h);

fork=1:j

D(j+1,k+1)=D(j+1,k)+(D(j+1,k)-D(j,k))/((4^k)-1);

15

end

err=abs(D(j+1,j+1)-D(j,j));

relerr=2*err/(abs(D(j+1,j+1))+abs(D(j,j))+eps);

j=j+1;

end

[n,n]=size(D);

【例】设)cos()(xxf，求)8.0('f的近似值，精度为小数点后9位。

f=inline('cos(x)');x=0.8;delta=0.5*10^-9;toler=delta;

[D,err,relerr,n]=diffext(f,x,delta,toler)

D=vpa(D,10),D_exact=vpa(-sin(0.8),10)

3、牛顿多项式微分（基于1n个点的差分求解）：通过构造下列n次牛顿多项式求解)('xf的近

似值：)()())(()()(

10102010



nn

xxxxaxxxxaxxaaxP

将)(')('

00

xPxf作为最终结果。在

0

x处使用这个方法。

可通过重新排列点的顺序为

nkkk

xxxxx,,,,,,

110





来计算)(')('

kk

xPxf。

function[A,df]=diffnew(X,Y)

%Input-Xisthe1xnabscissasvector

%-Yisthe1xnordinatevector

%Output-Aisthe1xnvectorcontainingthecoefficientsofthe

%Nth-degreeNewtonpolynomial

%-dfistheapproximatederivative

A=Y;

n=length(X);

forj=2:n

fork=n:-1:j

A(k)=(A(k)-A(k-1))/(X(k)-X(k-j+1));

end

end

x0=X(1);

df=A(2);

prod=1;

n1=length(A)-1;

fork=2:n1

prod=prod*(x0-X(k));

df=df+prod*A(k+1);

end

【例】X=0.8:0.05:1.0;Y=cos(X);[A,df]=diffnew(X,Y)

df=vpa(df,10),df_exact=vpa(-sin(0.8),10)

数值积分

1、复化梯形公式：通过)(xf的1n个等步长采样点：nkkhax

k

,,1,0,，其中

bxax

n

,

0

，逼近积分的数值公式为：







1

1

)()]()([

2

)(

n

k

k

b

a

xfhbfaf

h

dxxf。

functions=trap_rule(f,a,b,n)

%Input-fistheintegrandinputasastring'f'

%-aandbarelowerandupperlimitsofintegration

%-nisthenumberofsubintervals

16

%Output-sisthetrapezoidalrulesum

h=(b-a)/n;s=0;

fork=1:(n-1)

x=a+h*k;

s=s+feval(f,x);

end

s=h*(feval(f,a)+feval(f,b))/2+h*s;

【例】计算积分1

0

dxex。f=inline('exp(x)');a=0;b=1;n=5;

s=trap_rule(f,a,b,n);s=vpa(s,10)

积分准确值：exact=int('exp(x)',0,1);exact=vpa(exact,10)

2、复化辛普生公式：通过)(xf的12n个等步长采样点：nkkhax

k

2,,1,0,，其中

bxax

n



20

,，逼近积分的数值公式为：















n

k

k

n

k

k

b

a

xf

h

xf

h

bfaf

h

dxxf

1

12

1

1

2

)(

3

4

)(

3

2

)]()([

3

)(。

functions=trap_rule(f,a,b,n)

%Input-fistheintegrandinputasastring'f'

%-aandbarelowerandupperlimitsofintegration

%-nisthenumberofsubintervals

%Output-sisthetrapezoidalrulesum

h=(b-a)/n;s=0;

fork=1:(n-1)

x=a+h*k;

s=s+feval(f,x);

end

s=h*(feval(f,a)+feval(f,b))/2+h*s;

【例】计算积分1

0

dxex。a=0;b=1;n=1;s=simp_rule(f,a,b,n);s=vpa(s,10)

积分准确值：exact=int('exp(x)',0,1);exact=vpa(exact,10)

3、变步长梯形公式：由))()((

2

)0(bfaf

h

T开始，由如下递归公式可生成一个梯形公式

)}({jT的序列：,2,1,)(

2

)1(

)(

1

12











jxfh

jT

jT

n

k

k

，其中jabh2/)(，且

khax

k

。

functionT=rctrap(f,a,b,n)

%Input-fistheintegrandinputasastring'f'

%-aandbarelowerandupperlimitsofintegration

%-nisthenumberoftimesforrecursion

%Output-sistherecursivetrapezoidalrulelist

m=1;h=b-a;T=zeros(1,n+1);

T(1)=h*(feval(f,a)+feval(f,b))/2;

forj=1:n

m=2*m;

h=h/2;

s=0;

17

fork=1:m/2

x=a+h*(2*k-1);

s=s+feval(f,x);

end

T(j+1)=T(j)/2+h*s;

end

【例】计算积分5

1

1

dx

x

。

a=1;b=5;n=3;T=rctrap(f,a,b,n);T=vpa(T,10)

积分准确值：I_exact=int('1/x',1,5);I_exact=vpa(I_exact,10)

4、龙贝格积分：通过生成kj的逼近表),(kjR，并以)1,1(jjR为最终解来逼近积

分：

),()(jjRdxxfb

a

。

逼近),(kjR保存在一个特别的下三角矩阵中，第0列的元素)0,(jR用基于j2个],[ba子

区间的连续变步长梯形公式计算，然后利用龙贝格公式计算),(kjR。

第j行的元素为：jk

kjRkjR

kjRkjR

k







1,

14

)1,1()1,(

)1,(),(。

function[R,quad,err,h]=romberg(f,a,b,n,tol)

%Input-fistheintegrandinputasastring'f'

%-aandbarelowerandupperlimitsofintegration

%-nisthenumberofrowsinthetable

%-tolisthetolerance

%Output-RistheRombergtable

%-quadisthequadraturevalue

%-erristheerrorestimate

%-histhesmalleststepsizeud

m=1;h=b-a;err=1;j=0;R=zeros(4,4);

R(1,1)=h*(feval(f,a)+feval(f,b))/2;

while((err>tol)&(j

j=j+1;

h=h/2;

s=0;

fork=1:m

x=a+h*(2*k-1);

s=s+feval(f,x);

end

R(j+1,1)=R(j,1)/2+h*s;

m=2*m;

fork=1:j

R(j+1,k+1)=R(j+1,k)+(R(j+1,k)-R(j,k))/(4^k-1);

end

err=abs(R(j,j)-R(j+1,k+1));

end

quad=R(j+1,j+1);

【例】用龙贝格积分计算2/

0

2)cos()1(dxxxx

。

f=inline('(x^2+x+1)*cos(x)');a=0;b=pi/2;n=5;tol=0.0001;

[R,quad,err,h]=romberg(f,a,b,n,tol)

积分准确值：

I_exact=int('(x^2+x+1)*cos(x)',0,pi/2);I_exact=vpa(I_exact,10)

18

5、自适应辛普生积分

辛普生公式：子程序srule用于实现自适应积分过程中产生的子区间上的辛普生公式

))0()0(4)0((

3

)(0

0

bfcfaf

h

dxxfb

a

，其中2/)00(0bac。

functionZ=srule(f,a0,b0,tol0)

%Input-fistheintegrandinputasastring'f'

%-a0andb0arelowerandupperlimitsofintegration

%-tol0isthetolerance

%Output-Zisa1x6vector[a0b0ss2errtol1]

h=(b0-a0)/2;

c=zeros(1,3);

c=feval(f,[a0(a0+b0)/2b0]);

s=h*(c(1)+4*c(2)+c(3))/3;

s2=s;

tol1=tol0;

err=tol0;

Z=[a0b0ss2errtol1];

自适应辛普生积分：将对分前子区间],[

kk

ba细分为两个相等子区间],[

11kk

ba和],[

22kk

ba。

)]()(4)([

6

)]()(4)([

6

),(),(

2221112211kkkkkkkkkk

bfcfaf

h

bfcfaf

h

basbas

function[SRmat,quad,err]=adapt(f,a,b,tol)

%Input-fistheintegrandinputasastring'f'

%-aandbarelowerandupperlimitsofintegration

%-tolisthetolerance

%Output-SRmatisthetableofvalues

%-quadisthequadraturevalue

%-erristheerrorestimate

%Initialvalues

SRmat=zeros(30,6);

iterating=0;

done=1;

SRvec=zeros(1,6);

SRvec=srule(f,a,b,tol);

SRmat(1,1:6)=SRvec;

m=1;

state=iterating;

while(state==iterating)

n=m;

forj=n:-1:1

p=j;

SR0vec=SRmat(p,:);

err=SR0vec(5);

tol=SR0vec(6);

if(tol<=err)

%Bictinterval,applySimpson'srule

%recursively,anddetermineerror

state=done;

SR1vec=SR0vec;

SR2vec=SR0vec;

a=SR0vec(1);

b=SR0vec(2);

19

c=(a+b)/2;

err=SR0vec(5);

tol=SR0vec(6);

tol2=tol/2;

SR1vec=srule(f,a,c,tol2);

SR2vec=srule(f,c,b,tol2);

err=abs(SR0vec(3)-SR1vec(3)-SR2vec(3))/10;

%Accuracytest

if(err

SRmat(p,:)=SR0vec;

SRmat(p,4)=SR1vec(3)+SR2vec(3);

SRmat(p,5)=err;

el

SRmat(p+1:m+1,:)=SRmat(p:m,:);

m=m+1;

SRmat(p,:)=SR1vec;

SRmat(p+1,:)=SR2vec;

state=iterating;

end

end

end

end

quad=sum(SRmat(:,4));

err=sum(abs(SRmat(:,5)));

SRmat=SRmat(1:m,1:6);

【例】计算4

0

2/32)(13dxexxx。起始容差为001.0

0

。

该方法的实现需要6个子区间，输出项的6列依次为：每个子区间

kk

ba,，复化辛普生公

式对分前),(

kk

bas，对分后),(

kk

bas，该逼近的误差界以及相应的容差

k

。

f=inline('13.*(x-x.^2).*exp(-3*x./2)');a=0;b=4;tol=0.001;

[SRmat,quad,err]=adapt(f,a,b,tol),quad=vpa(quad,10)

积分准确值：

I_exact=int('13.*(x-x.^2).*exp(-3*x./2)',0,4);I_exact=vpa(I_exact,10)

6、高斯－勒让德积分:













n

k

kk

b

a

tfAdttf

ab

dxxf

1

1

1

)()(

2

)(,

其中

kk

x

abba

tx

abba

t

22

,

22















。

functionquad=guass(f,a,b,xk,Ak)

%Input-fistheintegrandinputasastring'f'

%-aandbarelowerandupperlimitsofintegration

%-xkisthe1xnvectorofabscissasfromTable

%-Akisthe1xnvectorofweightsfromTable

%Output-quadisthequadraturevalue

n=length(xk);

T=zeros(1,n);

T=(a+b)/2+((b-a)/2)*xk;

quad=((b-a)/2)*sum(Ak.*feval(f,T));

【例】利用两点高斯－勒让德积分计算







1

1

09861.1)1ln()3ln(

2

1

dx

x

。

f=inline('1./(x+2)');a=-1;b=1;xk=[-sqrt(3)/3,sqrt(3)/3];Ak=[11];

20

quad=guass(f,a,b,xk,Ak);quad=vpa(quad,10)

积分准确值：

I_exact=int('1/(x+2)',-1,1);I_exact=vpa(I_exact,10)

常微分方程数值方法

1、欧拉方法：

functionE=euler(f,a,b,ya,n)

%Input-fisthefunctionenteredasastring'f'

%-aandbaretheleftandrightendpoints

%-yaistheinitialconditiony(a)

%-nisthenumberofsteps

%Output-E=[T'Y']whereTisthevectorofabscissasand

%Yisthevectorofordinates

h=(b-a)/n;

T=zeros(1,n+1);

Y=zeros(1,n+1);

T=a:h:b;

Y(1)=ya;

forj=1:n

Y(j+1)=Y(j)+h*feval(f,T(j),Y(j));

end

E=[T'Y'];

【例】用欧拉方法求解区间]3,0[内的初值问题：1)0(,

2

'



y

yt

y。