《应用计算方法教程》matlab作业二Word格式.docx
《《应用计算方法教程》matlab作业二Word格式.docx》由会员分享，可在线阅读，更多相关《《应用计算方法教程》matlab作业二Word格式.docx（15页珍藏版）》请在冰豆网上搜索。
，令mk
max（Vk）,Uk
Vk
Vk1AUk,k1，当k
充分大时
mk
Uk
max（Vk1）
1。
max
（1）
QR法求全部特色值：
A
A1
Q1R1
RQ1
A2
Q2R2
k
1,2,3,
RkQk
Ak1
Qk1Rk1
由于此题的矩阵是10阶的，上述算法计算时间过长，考虑采用改良算法——移位加速。
迭代格式以下：
Ak
qkIQkRk
1RkQk
qkI
an（k
1,n）
an（k1,）
n
的特色值
（k）
（k）
，当
为实数时，选
qk为
计算Ak右下角的二阶矩阵
a
n1
n,n
n,n
n（k
1）,n（k）中最凑近an（k,n）的。
程序
A=-5+round（10*rand（10））;
[V,D]=eig（A）
[lamdau]=lab6_2_power（A,[1;
1;
1],10^（-5）,1000）
d=lab6_3_qr2（A,10^（-5））
function
lamda=0;
err=1;
[lamdau]=lab6_2_power（a,v,eps,N）
k=1;
while（k<
=N&
&
err>
eps）
u=a*v;
[mj]=max（abs（u））;
dc=abs（lamda-m）;
u=u/m;
dv=norm（u-v）;
err=max（dc,dv）;
v=u;
lamda=m;
k=k+1;
end
functionD=lab6_3_qr2（A,eps）
[n,n]=size（A）;
m=n;
D=zeros（n,1）;
B=A;
while（m>
1）
while（abs（B（m,m-1））>
=eps*（abs（B（m-1,m-1））+abs（B（m,m））））
S=eig（B（m-1:
m,m-1:
m））;
[j,k]=min（[abs（B（m,m）-S
（1））,abs（B（m,m）-S
（2））]）;
[Q,U]=qr（B-S（k）*eye（m））;
B=U*Q+S（k）*eye（m）;
A（1:
m,1:
m）=B;
m=m-1;
B=A（1:
m）;
D=diag（A）;
界面
（1）
（2）
（3）
作业七
7-1试验目的：
熟悉代数插值
在f（x）在7个点的函数值以下表所示，分别使用拉格朗日插值法和牛顿插值法求f（0.596）与f（0.906）的近似值。
xi
yi
拉格朗日插值多项式：
（xx0）（xx1）（xxj1）（xxj1）
（xxn）
Ln（x）
j
（xj
x0）（xjx1）（xjxj1）（xjxj1）
（xjxn）yj
（
x
xi）yj
i
xj
牛顿插值多项式：
Nn（x）
f[x0,
f（x0）,xn]（x
f[x0,x1]（x
x0）（x
x0xn1
）
f[x0,x1,x2
]（x
x0）（x
x1）
f（xj
其中f[x0,,xk]
。
j0（xjx0）
（xjxj1）（xjxj1）（xjxk）
functiony1=lab7_1_Lagrange（x,y,x1）
y1=0;
[mn]=size（x）;
n=n-1;
fork=1:
n+1
t=1;
fori=1:
if（i~=k）
t=t*（x1-x（i））/（x（k）-x（i））;
y1=y1+t*y（k）;
functiony1=lab7_2_Newton（x,y,x1）
forj=1:
fori=n+1:
-1:
j+1
y（i）=（y（i）-y（i-1））/（x（i）-x（i-j））;
y1=y（n+1）;
forj=n:
y1=y（j）+（x1-x（j））*y1;
作业八
8-1试验目的：
熟悉最小二乘法拟合多项式
给定数据点（xi,yi），
X
f（x）
用3次最小二乘多项式拟合数据，并求平方误差。
要作三次最小二乘拟合，令
S（x）
a0
0（x）
a11（x）a22（x）
a3
3（x）,
j（x）xj，计算法方程
（0,0）
（0,1）
（0,n）
d1
Gad，其中G
1,0）
1,
1）
1,n），a
a1
，d
d2
，
（n,0）
（n,1）
（n,n）
an
dn
m
j,k）
（xi
j（xi）
k（xi）,（j,k
0,1,,n）
（xi）
dk
（f,k）
（xi）f（xi）k（xi）,（k
0,1,
n）
其平方误差为s
f（xi）]2。
i0
（xi）[S（xi
x=[0.40.550.650.800.901.05];
f=[0.410750.578150.696750.888111.026521.25386];
G=zeros（4,4）;
4
6
G（j,k）=G（j,k）+x（i）^（j+k-2）;
d=zeros（4,1）;
d（k,1）=d（k,1）+f（i）*x（i）^（k-1）;
a=G\d
s=0;
s=s+（f（i）-（a
（1）+a
（2）*x（i）+a（3）*x（i）^2+a（4）*x（i）^3））^2;
s
作业九
9-1试验目的：
熟悉数值积分公式，掌握数值计算定积分的方法
采用不同样方法数值计算积分
1ln（1x）
dx
编写复合梯形公式和复合
Simpson公式通用子程序，分别采用
4，8，16，32，64均分区间
计算。
复合梯形公式：
将区间
[a,b]作n均分，h
ba，结点xi
ih（0in），
Tn
n1h
h
f（xi）
f（b）]
[f（xi）f（xi1）]
[f（a）
2
i02
i1
复合Simpson公式：
将区间[a,b]作2n均分，记h
b
a，
2n
S2n
h[f（x2i）4f（x2i1）f（x2i2）]
h[f（a）
f（x2i1）
f（x2i）f（b）]
3
functiony=lab9_f（x）
y=（log（1+x））/x;
functiony=lab9_1_fTrapezoid（a,b,eps,n）
f0=0;
h=（b-a）/n;
fori=0:
n-1
f0=f0+lab9_f（a+i*（b-a）/n）+lab9_f（a+（i+1）*（b-a）/n）;
y=f0*h/2;
functiony=lab9_2_fSimpson（a,b,eps,n）
k=n/2;
k-1
f0=f0+lab9_f（a+2*i*（b-a）/n）+4*lab9_f（a+（2*i+1）*（b-a）/n）+lab9_f（a+（2*i+2）*（b-a）/n）;
y=f0*h/3;
作业十
10-1试验目的：
学会用Euler法、改良Euler法、经典的4阶Runge-Kutta法求解常微分方程
初值问题。
分别用
1）Euler法〔步长〕
2）改良Euler法〔步长〕
3）4阶Runge-Kutta〔步长〕
求解下面的初值问题：
y
（y
1）（y
3）,0x
比较公共节点解的误差。
精确解为
2（1
e2x）1。
Euler法：
令f（x,y（x））
y（x），yn
yn
hf（xn
yn）〔n
ba〕。
k1
f（xn,yn）
改良Euler
法〔梯形公式〕：
k2
f（xn
h,yn
hk1）
h（k
hf（xn,yn）
k2
hf（xn
yn
k1）
4阶Runge-Kutta：
k3
h,yn
1k2）
k4
k3）
yn1
2k2
2k3
k4）
（k1
y0=-2;
h=0.025;
n=2/h;
y
（1）=y0-h*（y0+1）*（y0+3）;
fori=2:
y（i）=y（i-1）-h*（y（i-1）+1）*（y（i-1）+3）;
y1（j）=-3+2/（1+exp（-4*j/n））;
y1
x=0.025:
0.025:
2;
plot（x,y1,'
r'
x,y,
'
b'
h=0.05;
k1=-（y0+1）*（y0+3）;
k2=-（y0+h*k1+1）*（y0+h*k1+3）;
y21
（1）=y0+h/2*（k1+k2）;
k1=-（y21（i-1）+1）*（y21（i-1）+3）;
k2=-（y21（i-1）+h*k1+1）*（y21（i-1）+h*k1+3）;
y21（i）=y21（i-1）+h/2*（k1+k2）;
y21
y22（j）=-3+2/（1+exp（-4*j/n））;
y22
x=0.05:
0.05:
plot（x,y21,
x,y22,
h=0.1;
k1=-h*（y0+1）*（y0+3）;
k2=-h*（y0+k1/2+1）*（y0+k1/2+3）;
k3=-h*（y0+k2/2+1）*（y0+k2/2+3）;
k4=-h*（y0+k3+1）*（y0+k3+3）;
y31
（1）=y0+（k1+2*k2+2*k3+k4）/6;
k1=-h*（y31（i-1）+1）*（y31（i-1）+3）;
k2=-h*（y31（i-1）+k1/2+1）*（y31（i-1）+k1/2+3）;
k3=-h*（y31（i-1）+k2/2+1）*（y31（i-1）+k2/2+3）;
k4=-h*（y31（i-1）+k3+1）*（y31（i-1）+k3+3）;
y31（i）=y31（i-1）+（k1+2*k2+2*k3+k4）/6;
y31
y32（j）=-3+2/（1+exp（-4*j/n））;
y32
x=0.1:
0.1:
plot（x,y31,
x,y32,
升级会员 