矩阵与数值分析Word下载.docx
《矩阵与数值分析Word下载.docx》由会员分享，可在线阅读，更多相关《矩阵与数值分析Word下载.docx（25页珍藏版）》请在冰点文库上搜索。
T（i,j）=2;
end
if（i==（j+1））
T（i,j）=-1;
if（i==（j-1））
I=eye（n-1）;
k=1;
formm=1:
（n-1）
A（k:
（k+n-2）,k:
（k+n-2））=T+2*I;
k=k+n-1;
forxx=1:
（n-2）
（k+n-2）,（k+n-1）:
（k+2*n-3））=-I;
A（（k+n-1）:
（k+2*n-3）,k:
（k+n-2））=-I;
b=randn（（n-1）^2,1）;
x0=zeros（（n-1）^2,1）;
Jacobi迭代法Matlab程序如下：
functionn=jacobi（A,b,x0）
eps=1.0e-10;
M=10000;
D=diag（diag（A））;
%求A的对角矩阵
L=-tril（A,-1）;
%求A的下三角阵
U=-triu（A,1）;
%求A的上三角阵
B=D\（L+U）;
f=D\b;
x=B*x0+f;
n=1;
%迭代次数
whilenorm（x-x0）>
=eps
x0=x;
x=B*x0+f;
n=n+1;
if（n>
=M）
disp（'
Warning:
迭代次数太多，可能不收敛！
'
）;
return;
Gauss-Seidel迭代法Matlab程序如下：
functionn=gauseidel（A,b,x0）
G=（D-L）\U;
f=（D-L）\b;
x=G*x0+f;
x=G*x0+f;
[A,b,x0]=jz（10）;
J10=jacobi（A,b,x0）
G10=gauseidel（A,b,x0）
[A,b,x0]=jz（20）;
J20=jacobi（A,b,x0）
G20=gauseidel（A,b,x0）
[A,b,x0]=jz（30）;
J30=jacobi（A,b,x0）
G30=gauseidel（A,b,x0）
J10=436；
G10=214
J20=1810；
G20=913
J30=3990；
G30=2001
N=10且M=1000时，Jacobi迭代法和Gauss—seidel迭代法的迭代次数分别为436和214；
N=20且M=1000时，Jacobi迭代法和Gauss—seidel迭代法的迭代次数分别为1810和913；
N=30且M=1000时，Jacobi和Gauss-seidel算法的迭代次数分别为3990和2001次。
从以上结果可知在该问题下Jacobi算法的收敛性没有Gauss-seidel算法的收敛性好，原因在于Jacobi算法迭代过程中，对已算出的信息未加充分利用，一般来说，后面计算的值要比前面的计算值要精确些，而Gauss-seidel算法则充分利用了已求出来的信息，所以此算法的收敛性更为好一些。
2.用Gauss列主元消去法、QR方法求解如下方程组：
Gauss列主元消去法Matlab程序如下：
function[x,XA]=GaussXQLineMain（A,b）
N=size（A）;
n=N
（1）;
index=0;
me=max（abs（A（1:
n,i）））;
%选取列主元
fork=i:
n
if（abs（A（k,i））==me）
index=k;
%保存列主元所在的行
break;
temp=A（i,1:
n）;
A（i,1:
n）=A（index,1:
A（index,1:
n）=temp;
bb=b（index）;
b（index）=b（i）;
b（i）=bb;
%交换主行
forj=（i+1）:
if（A（i,i）==0）
对角元素为0！
l=A（j,i）;
m=A（i,i）;
A（j,1:
n）=A（j,1:
n）-l*A（i,1:
n）/m;
b（j）=b（j）-l*b（i）/m;
%消元
fori=n:
1
if（i<
n）
s=A（i,（i+1）:
n）*x（（i+1）:
n,1）;
else
s=0;
x（i,1）=（b（i）-s）/A（i,i）;
end
XA=A;
QR法Matlab程序如下：
functionx=qrxq（A,b）
B=A;
%保存系数矩阵
A（1:
n,1）=A（1:
n,1）/norm（A（1:
n,1））;
%将A的第一列正规化
fori=2:
（i-1）
A（1:
n,i）=A（1:
n,i）-dot（A（1:
n,i）,A（1:
n,j））*A（1:
n,j）;
%使A的第i列与前面所有的列正交
n,i）=A（1:
n,i）/norm（A（1:
n,i））;
%将A的第i列正规化
Q=A;
%分解出来的正交矩阵
R=transpose（Q）*B;
bb=transpose（Q）*b;
s=R（i,（i+1）:
x（i,1）=（bb（i）-s）/R（i,i）;
A=[12-11;
2503;
1792;
8-1-21];
b=[0;
4;
12;
-8];
x_G=GaussXQLineMain（A,b）
x_QR=qrxq（A,b）
x_G=-1.00000000000000
-0.00000000000000
1.00000000000000
2.00000000000000
x_QR=-1.00000000000000
-0.00000000000001
1.00000000000000
2.00000000000002
三、非线性方程的迭代解法
1．求方程
的根，迭代停止的条件为：
；
使用弦截法求解此方程的根，弦截法的Matlab程序如下：
functionroot=Secant（f,a,b）
%f:
方程;
%a:
区间左端点;
%b:
区间右端点;
%root:
求得的根;
eps=1.0e-10;
f1=subs（sym（f）,findsym（sym（f））,a）;
f2=subs（sym（f）,findsym（sym（f））,b）;
if（f1==0）
root=a;
if（f2==0）
root=b;
if（f1*f2>
0）
两端点函数值乘积大于0!
else
tol=1;
fa=subs（sym（f）,findsym（sym（f））,a）;
fb=subs（sym（f）,findsym（sym（f））,b）;
root=a-（b-a）*fa/（fb-fa）;
while（tol>
eps）
r1=root;
fx=subs（sym（f）,findsym（sym（f））,r1）;
s=fx*fa;
if（s==0）
root=r1;
if（s>
root=b-（r1-b）*fb/（fx-fb）;
root=a-（r1-a）*fa/（fx-fa）;
tol=abs（root-r1）;
利用matlab的绘图功能可以确定求根区间为（1,3）
root=Secant（'
exp（x）+2*x^2+2*sin（x）-log（x）-16'
1,3）
root=1.96292268316446
2．利用Newton迭代法求多项式
的所有实零点，注意重根的问题。
Newton迭代法的Matlab程序如下：
x=-3;
d=-1;
whileabs（d-x）>
0.5*10^（-8）%在区间（-3，-1）
d=x;
x=x-（x^4-3*x^3-3*x^2+11*x-6）/（4*x^3-9*x^2-6*x+11）;
%牛顿法
x0=d
x=0;
d=2;
0.5*10^（-8）%在区间（0，2）的重根
x=x-2*（x^4-3*x^3-3*x^2+11*x-6）/（4*x^3-9*x^2-6*x+11）;
x1=d
x=2.5;
d=4;
0.5*10^（-8）%在区间（2.5，4）
x2=d
x0=-2.00000000151233
x1=1.00000000065769
x2=3.00000000000002
由以上结果可知，newton迭代法所得结果与初始值的选取密切相关。
不同的初始值会产生不同的结果，甚至会影响牛顿迭代法收敛性。
四、数值积分
用三点Gauss型求积公式计算积分
n点Gauss型求积公式matlab程序如下：
functionI=question4（f,a,b,n）
%f为求积方程;
%a,b分别为积分上下限;
%n为Gauss点个数;
[XK,AK]=grule（n）;
XK1=（（b-a）/2）*XK+（（a+b）/2）;
I=（（b-a）/2）*sum（AK.*subs（f,findsym（f）,XK1））;
function[bp,wf]=grule（n）
%[bp,wf]=grule（n）
bp=zeros（n,1）;
wf=bp;
iter=2;
m=fix（（n+1）/2）;
e1=n*（n+1）;
mm=4*m-1;
t=（pi/（4*n+2））*（3:
4:
mm）;
nn=（1-（1-1/n）/（8*n*n））;
xo=nn*cos（t）;
forj=1:
iter
pkm1=1;
pk=xo;
fork=2:
t1=xo.*pk;
pkp1=t1-pkm1-（t1-pkm1）/k+t1;
pkm1=pk;
pk=pkp1;
den=1.-xo.*xo;
d1=n*（pkm1-xo.*pk）;
dpn=d1./den;
d2pn=（2.*xo.*dpn-e1.*pk）./den;
d3pn=（4*xo.*d2pn+（2-e1）.*dpn）./den;
d4pn=（6*xo.*d3pn+（6-e1）.*d2pn）./den;
u=pk./dpn;
v=d2pn./dpn;
h=-u.*（1+（.5*u）.*（v+u.*（v.*v-u.*d3pn./（3*dpn））））;
p=pk+h.*（dpn+（.5*h）.*（d2pn+（h/3）.*（d3pn+.25*h.*d4pn）））;
dp=dpn+h.*（d2pn+（.5*h）.*（d3pn+h.*d4pn/3））;
h=h-p./dp;
xo=xo+h;
bp=-xo-h;
fx=d1-h.*e1.*（pk+（h/2）.*（dpn+（h/3）.*（...
d2pn+（h/4）.*（d3pn+（.2*h）.*d4pn））））;
wf=2*（1-bp.^2）./（fx.*fx）;
if（m+m）>
n,bp（m）=0;
if~（（m+m）==n）,m=m-1;
jj=1:
m;
n1j=（n+1-jj）;
bp（n1j）=-bp（jj）;
wf（n1j）=wf（jj）;
%end
symsx;
f=exp（-x）;
a=0;
b=1;
n=3;
I=question4（f,a,b,n）
I=0.632120255664068
五、插值与逼近
1．给定
上的函数
，请做如下的插值逼近：
（1）构造等距节点分别取
的Lagrange插值多项式；
（2）取Chebyshev多项式
的零点：
作插值节点构造
的插值多项式
和上述的插值多项式均要求画出曲线图形（用不同的线型或颜色表示不同的曲线）。
（1）Lagrange插值matlab程序如下：
functionf=L1（n）
symsts;
q=1/（1+s^2）;
x（i）=-1+2*（i-1）/（n-1）;
y（i）=subs（q,s,x（i））;
f=0.0;
for（i=1:
l=y（i）;
for（j=1:
i-1）
l=l*（t-x（j））/（x（i）-x（j））;
end;
for（j=i+1:
%计算拉格朗日基函数
f=f+l;
%计算拉格朗日插值函数
f=collect（f）;
%化简
X=-1:
0.1:
1;
Y=subs（q,s,X）;
F=subs（f,t,X）;
plot（X,Y,'
g'
X,F,'
r'
f=collect（f）;
%展开
digits（5）;
f=vpa（f）;
1、在matlab命令窗口中输入：
f=L1（5）
运行结果
2、在matlab命令窗口中输入：
f=L1（8）
3、在matlab命令窗口中输入：
f=L1（10）
由图像可以看出n=5时拟合曲线与原始曲线的误差比较大，而n=8,10时拟合曲线与原始曲线几乎重合在一起。
可见节点个越多，拉格朗日插值法所得的多项式曲线与原始曲线越接近，误差就越小。
的零点为插值点,此时Lagrange插值matlab程序主体基本不变，变的只是插值点，其程序如下：
functionf=L2（n）
x（i）=cos（（2*i-1）*pi/2/n）;
%计算拉格朗日插值函数
%绘图
k'
rp'
f=L2（10）
有图像比较可以看出，N=10时，chebyshev零为拉格朗日插值时，其拟合曲线与原始曲线几乎重合，精度较高。
2．已知函数值
3
4
6
和边界条件：
． 求三次样条插值函数
并画出其图形．
三次样条插值matlab程序如下：
functionThrSample1（x,y,y_1,y_N）
symst;
if（length（x）==length（y））
n=length（x）;
x和y的维数不相等！
end%维数检查
A=diag（2*ones（1,n））;
%求解m的系数矩阵
u=zeros（n-2,1）;
lamda=zeros（n-1,1）;
c=zeros（n,1）;
u（i-1）=（x（i）-x（i-1））/（x（i+1）-x（i-1））;
lamda（i）=（x（i+1）-x（i））/（x（i+1）-x（i-1））;
c（i）=3*lamda（i）*（y（i）-y（i-1））/（x（i）-x（i-1））+...
3*u（i-1）*（y（i+1）-y（i））/（x（i+1）-x（i））;
A（i,i+1）=u（i-1）;
A（i,i-1）=lamda（i）;
%形成系数矩阵及向量c
c
（1）=2*y_1;
c（n）=2*y_N;
m=followup（A,c）;
%用追赶法求解方程组
j=1;
whilej<
=n-1
h=x（j+1）-x（j）;
f=y（j）*（2*（t-x（j））+h）*（t-x（j+1））^2/h^3+...
y（j+1）*（2*（x（j+1）-t）+h）*（t-x（j））^2/h^3+...
m（j）*（t-x（j））*（x（j+1）-t）^2/h^2-...
m（j+1）*（x（j+1）-t）*（t-x（j））^2/h^2;
fprintf（'
当x属于[%d,%d],f=\n'
x（j）,x（j+1））;
disp（f）;
X=x（j）:
x（j+1）;
F=subs（f,t,X）;
plot（X,F,'
holdon;
j=j+1;
functionx=followup（A,b）
n=rank（A）;
for（i=1:
Error:
对角有元素为0！
end;
d=ones（n,1）;
a=ones（n-1,1）;
c=ones（n-1）;
n-1）
a（i,1）=A（i+1,i）;
c（i,1）=A（i,i+1）;
d（i,1）=A（i,i）;
d（n,1）=A（n,n）;
for（i=2:
n）
d（i,1）=d（i,1）-（a（i-1,1）/d（i-1,1））*c（i-1,1）;
b（i,1）=b（i,1）-（a（i-1,1）/d（i-1,1））*b（i-1,1）;
x（n,1）=b（n,1）/d（n,1）;
for（i=（n-1）:
1）
x（i,1）=（b（i,1）-c（i,1）*x（i+1,1））/d（i,1）;
x=[346];
y=[602];
y_1=1;
y_N=-1;
ThrSample1（x,y,y_1,y_N）
由图像可以看出三次样条插值曲线的优越性，即光滑性好，收敛得到保证，局部性好。
3．观察物体的直线运动，得到如下数据
时刻t
0.9
1.9
3.0
3.9
5.0
位移s
10
30
51
80
111
求运动学方程
。
matlab程序如下：
t=[00.91.93.03.95.0];
s=[010305180111];
n=length（t）;
a2=0;
a2=t（i）^2+a2;
a3=0;
a3=t（i）^3+a3;
a4=0;
a4=t（i）^4+a4;
b1=0;
b1=s（i）.*t（i）+b1;
b2=0;
b2=s（i）.*（t（i）.^2）+b2;
A=[a2a3;
a3a4]
b=[b1;
b2]
c=inv（A）*b
x=0:
0.01:
5.0;
y=c
（1）.*x+c（
升级会员 