常州大学数值分析课后习题答案第二章第三章第四章节资料文档格式.docx
《常州大学数值分析课后习题答案第二章第三章第四章节资料文档格式.docx》由会员分享，可在线阅读，更多相关《常州大学数值分析课后习题答案第二章第三章第四章节资料文档格式.docx（16页珍藏版）》请在冰点文库上搜索。
x
A=2-131
4254
1207
x=9
-1
-6
11x1-3x2-2x3=3,
（2）-23x1+11x2+1x3=0,
x1+2x2+2x3=-1;
（2）解：
A=[11-3-23;
-231110;
122-1]
%消元过程
n-1
n
eps
n+1）=
return
end
%回代过程
x（n）=A（n,n+1）/A（n,n）;
else
1
x（i）=（A（i,n+1）-A（i,i+1:
x
A=11-3-23
-231110
122-1
x=0.2124
0.5492
-1.1554
4、用Cholesky分解法解方程组
323x15
220x23
3012x37
解：
.
A=[323;
220;
3012];
b=[537];
lambda=eig（A）;
iflambda>
eps&
isequal（A,A'
）
[n,n]=size（A）;
R=chol（A）;
%解R'
y=b
y
（1）=b
（1）/R（1,1）;
ifn>
fori=2:
y（i）=（b（i）-R（1:
i-1,i）'
*y（1:
i-1）'
）/R（i,i）;
%解Rx=y
x（n）=y（n）/R（n,n）;
x（i）=（y（i）-R（i,i+1:
n）'
x=x'
;
x=[];
disp（'
该方法只适用于对称正定的系数矩阵!
'
）;
R=1.73211.15471.7321
00.8165-2.4495
001.7321
y=2.8868-0.40820.5774
x=1.00000.50000.3333
5.用列主元Doolittle分解法解方程组
A=[345;
-134;
-23-5;
];
345X12
b=[2,-26]'
-134X2-2
[L,U,pv]=luex（A）;
-23-5X36
y=L\b（pv）;
x=U\y
结果如下：
x= 1
1
-1
14.已知，计算
.
A=[10099;
9998];
cond（A,inf）
ans=3.9601e+04
cond（A,2）
ans=3.9206e+04
27.编写LU分解法，改进平方根法，追赶法的Matlab程序，并进行相关数值试验。
LU分解法程序
Function[L,U]=lup（A）
%lup:
LUfactorization
%Synopsis:
[L,U]=lup（A）
%Input:
A=coefficientmatrix
%Output:
L:
lowertriangularmatrix
%Uuppertriangularmatrix
Formatshort
[m,n]=size（A）;
Ifm~=n,error（`Amatrixneedstobesquare`）;
End
Pv=（1:
n）`;
%LUfactorization
Fori=1:
Pivot=A（i,i）;
Fork=i+1;
A（k,i）=A（k,i）/pivot;
A（k,i+1;
n）=A（k,i+1;
n）-A（k,i）*A（i,i+1;
n）;
L=eye（size（A））+tril（A,-1）;
%extractLandU
U=triu（A）
改进平方根法程序
Function[x]=ave（A,b,n）
L=zeros（n,n）;
D=diag（n,0）;
S=L*D;
L（i;
i）=1;
Fori=1;
Forj=1;
If（eig（A）<
=0）|（A（i,j）~=A（j,i））disp（`wrong`）;
Break;
D（1,1）=A（1,1）；
Fori=2;
i-1
S（i,j）=A（i,j）-sum（S（i,1;
i-1）*L（j,1;
j-1）`）;
L（i,1;
i-1）=S（i,1;
i-1）/D（1;
i-1,1;
i-1）;
D（i,i）=A（i,i）-sum（L（i,1;
i-1）*L（i,1;
i-1）`）;
D（i,i）=A（i,i）-sum（S（i,1;
i-1）*D（1;
i-1）*y（1;
i-1）））/D（i,i）;
Y=zeros（n,1）;
X=zero（n,1）;
Y（i）=（b（i）-sum（L（i,1;
Fori=n;
-1;
X（i）=y（i）-sum（L（i+1;
n,i）`*x（i+1;
n））;
追赶法程序
Function[x,L,U]=Thomas（a,b,c,f）
N=length（b）;
%对A进行分解
U
（1）=b
（1）;
If（u（i-1）`=0）
L（i-1）=a（i-1）/u（i-1）;
U（i）=b（i）-l（i-1）*c（i-1）;
Else
L=eye（n）+diag（1,-1）;
U=diag（u）+diag（c,1）;
X=zeros（n,1）;
Y=x;
%?
求解ly=b?
Y
（1）=f
（1）;
Y（i）=f（i）-l（i-1）*y（i-1）;
求解Ux=y?
If（u（n）`=0）
X（n）=y（n）/u（n）;
Fori=n-1;
X（i）=（y（i）-c（i）*x（i+1））/u（i）；
第三章
1、设节点x0=0,x1=π/8,x2=π/4,x3=3π/8,x4=π/2,适当选取上述节点用Lagrange插值法分别构造cosx在区间[0,π/2]上的一次，二次和四次插值多项式P1（x）P2（x）和P4（x）,并分别计算P1（x）,P2（x）,P4（x）其中X取π/3。
A=fliplr（A）;
Return
x=[π/8,3π/8];
y=cos（x）;
x0=π/3;
[A,Y]=lagrange（x,y,x0）;
P1=vpa（poly2sym（A）,3）Y
P1=1.19x-0.689Y=0.4729x0=π/3;
P2=vpa（poly2sym（A）,3）Y
P2=x2-0.109x-0.336Y=0.5174
x=[0,π/8,π/4,3π/8,π/2];
y=cos（x）;
[A,Y]=lagrange（x,y,x0）;
P4=vpa（poly2sym（A）,3）Y
P4=x4+0.00282x3-0.514x2+0.0232x+0.0287Y=0.5001
7.根据列表函数
0.0
1.0
2.0
3.0
4.0
f（x）
0.50
1.25
2.75
3.50
选取适当的节点，用逐次线性插值法给出三次多项式在2.8处的值。
答：
Matlab程序function
[T,y0]=aitken（x,y,x0,T0）ifnargin==3T0=[];
end
n0=size（T0,1）;
m=max（size（x））;
n=n0+m;
T=zeros（n,n+1）;
T（1:
n0,1:
n0+1）=T0;
T（n0+1:
n,1）=x;
T（n0+1:
n,2）=y;
ifn0==0i0=2;
else
i0=n0+1;
fori=i0:
forj=3:
i+1
T（i,j）=fun（T（j-2,1）,T（i,1）,T（j-2,j-1）,T（i,j-1）,x0）;
end
y0=T（n,n+1）;
return
function[y]=fun（x1,x2,y1,y2,x）y=y1+（y2-y1）*（x-x1）/（x2-x1）;
%选取0、1、3、4四个节点，求三次插值多项式x=[0,1,3,4];
y=[0.5,1.25,3.5,2.75];
x0=2.8;
[T,y0]=aitken（x,y,x0）
y0=3.419000000000000
8.根据上题中的列表函数，写出差商表，并写出Newton插值多项式N2（x）和N4（x）。
差商表：
0.5
0.75
1.125
0.375
0.125
-0.25
0.5625
-0.0625
-0.21875
0.03125
由Nn（x）=a0+a1（x-x0）+a2（x-x0）（x-x1）+…+an（x-x0）（x-x1）…（x-xn-1）得
N2（x）=0.5+0.75（x-0.00）+0.375（x-0.0）（x-1.0）
=0.375x2+0.375x+0.5
N4（x）=0.5+0.75（x-0.00）+0.375（x-0.0）（x-1.0）-0.25（x-0.0）（x-1.0）（x-2.0）+0.03125（x-0.0）（x-2.0）（x-1.0）（x-3.0）
=0.03125x4-0.4375x3+1.46875x2-0.3125x+0.5
16、选取适当的函数y=f（x）和插值节点，编写Matlab程序,分别利用Lagrange插值方法,Newton插值方法确定的插值多项式，并将函数y=f（x）的插值多项式和插值余项的图形画在同一坐标系中，观测节点变化对插值余项的影响。
[C,D,Y]=newpoly（x0,y0,x）
%检验输入参数
ifnargin<
2|nargin>
3error（'
IncorrectNumberofInputs'
iflength（x0）~=length（y0）error（'
Thelengthofx0mustbeequaltoitofy0'
n=length（x0）;
D=zeros（n,n）;
D（:
1）=y0'
%计算差商表
forj=2:
nfork=j:
ifabs（x0（k）-x0（k-j+1））<
epserror（'
DividedbyZero,therearetwonodesarethesame'
D（k,j）=（D（k,j-1）-D（k-1,j-1））/（x0（k）-x0（k-j+1））;
endend
%计算Newton插值多项式的系数C=D（n,n）;
fork=（n-1）:
C=conv（C,poly（x0（k）））;
m=length（C）;
C（m）=C（m）+D（k,k）;
ifnargin==3
Y=polyval（C,x）;
x=[01234];
y=[0.5,1.25,2.75,3.5,2.75];
x0=[01234];
y0=[0.5,1.25,2.75,3.5,2.75];
%用lagrang插值法计算[A,Y]=lagrange（x,y,x0）
Lx=vpa（poly2sym（A）,4）%用newton插值法计算[C,D,X]=newpoly（x0,y0,x）Nx=vpa（poly2sym（C）,4）%绘制两者图像plot（x,Y,'
b*'
x0,X,'
r-'
）
A=0.5000-0.31251.4687-0.43750.0313
Y=0.50001.25002.75003.50002.7500
Lx=0.5x4-0.3125x3+1.469x2-0.4375x+0.03125
C=0.5000-0.31251.4688-0.43750.0313
D=0.50000000
1.25000.7500000
2.75001.50000.375000
3.50000.7500-0.3750-0.25000
2.7500-0.7500-0.7500-0.12500.0313
X=0.50001.25002.75003.50002.7500
Nx=0.5x4-0.3125x3+1.469x2-0.4375x+0.0312
第四章
6、已知列表函数
x1.02.03.04.05.0
y1.2222.9845.4668.90213.592
用最小二乘法求形如y=axebx的拟合函数。
Matlab程序
function[a,b]=ec（x,y）Y=log（y）'
A=zeros（5,3）;
fori=1:
5A（i,1）=1;
A（i,2）=log（x（i））;
A（i,3）=i;
c=inv（A'
*A）*（A'
*Y）;
a=exp（c
（1））;
b=c（3）;
5
y=a*x.*exp（b*x）;
endreturn
x=[12345];
y=[1.2222.9845.4668.90213.592];
[a,b]=ec（x,y）输出结果为：
a=1.000202219673205b=0.200293860504786
8、学习Matlab内部的函数lsqcurvefit，并设计数值实验使用lsqcurvefit。
Matlab内部函数lsqcurvefit是用来解决非线性拟合的最小二乘问题的。
其调用格式为：
x=lsqcurvefit（fun,x0,xdata,ydata）x=lsqcurvefit（fun,x0,xdata,ydata,lb,ub）
x=lsqcurvefit（fun,x0,xdata,ydata,lb,ub,options）
[x,resnorm]=lsqcurvefit（…）
[x,resnorm,residual]=lsqcurvefit（…）
[x,resnorm,residual,exitflag]=lsqcurvefit（…）
[x,resnorm,residual,exitflag,output]=lsqcurvefit（…）
[x,resnorm,residual,exitflag,output,lambda]=lsqcurvefit（…）[x,resnorm,residual,exitflag,output,lambda,jacobian]=lsqcurvefit（…）
输入参数：
fun为待拟合函数，计算x处拟合函数值，其定义为functionF=myfun（x,xdata）
x0为初始解向量，即拟合参数的初始解；
xdata，ydata为满足关系ydata=F（x,xdata）的数据；
lb、ub为解向量的下界和上界lb≤x≤ub，若没有指定界，则lb=[]，ub=[]；
options为指定的优化参数；
输出参数：
x为迭代得出解向量，即拟合出的参数；
resnorm=sum（（fun（x,xdata）-ydata）.^2），即x处残差平方和,最小二乘式值;
residual=fun（x,xdata）-ydata，即在x处的残差；
exitflag为终止迭代的条件；
output为输出的优化信息；
lambda为解x处的Lagrange乘子；
jacobian为解x处拟合函数fun的jacobian矩阵。
functionF=myfun（x,xdata）
F=（x
（1）.*xdata）.*（exp（x
（2）.*xdata））;
xdata=[1,2,3,4,5];
ydata=[1.222,2.984,5.466,8.902,13.592];
x0=[0,0];
[x,resnorm]=lsqcurvefit（@myfun,x0,xdata,ydata）
输出结果为：
Localminimumfound.
Optimizationcompletedbecausethesizeofthegradientislessthanthedefaultvalueofthefunctiontolerance.
<
stoppingcriteriadetails>
x=
0.9999583489763910.200014132812834resnorm=8.067930437509675e-7
升级会员 