北科计算方法上机作业Word下载.docx
《北科计算方法上机作业Word下载.docx》由会员分享，可在线阅读，更多相关《北科计算方法上机作业Word下载.docx（21页珍藏版）》请在冰豆网上搜索。
N
p1=p0-（3*p0-5^p0+4）/（3-5^p0*log（5））;
ifabs（p1-p0）<
Tol
break
p0=p1;
disp（p1）;
近似解以及迭代次数
当初值p0取-1时，迭代公式不变，仍是
pk+1=pk-[（3pk-5pk+4）/（3-5pk*ln5）]，程序不变（只是p0=-1）
结果如下：
结果分析：
当初值p0=1时，设不动点每步迭代的结果为p牛顿迭代法每步迭代的结果为p1，则
n
1
2
3
4
5
6
7
p
1.2091
1.2624
1.2753
1.2784
1.2791
1.2793
p1
1.3963
1.2922
1.2795
当初值p0=-1时，设不动点每步迭代的结果为p´
，牛顿迭代法每步迭代的结果为p1´
，则
p´
-1
-1.2667
-1.2899
-1.2915
-1.2916
p1´
-1.2987
由以上的迭代结果可以看出，Newton法迭代相对于不动点迭代法，收敛速度更快。
2、用Newton法与重根计算法求解方程x-sinx＝0的根.再用Steffensen’smethod加速Newton法收敛，比较结果。
A）Newton法求解：
迭代公式pk+1=pk-[pk-sin（pk）]/（1-cospk）
p1=p0-（p0-sin（p0））/（1-cos（p0））;
disp（k）近似解以及迭代次数：
B）重根法求解：
设f（x）=x-sinx，则有f（0）=f´
（0）=f´
´
（0），f（3）（0）=1，则p=0是方程的三重根。
于是有迭代公式pk+1=pk-3[pk-sin（pk）]/（1-cospk）
p1=p0-3*（p0-sin（p0））/（1-cos（p0））;
C）Steffensen’smethod加速Newton法收敛：
p1=0.6551;
n=0;
p
（1）=p0;
whilen<
p
（2）=p
（1）-（p
（1）-sin（p
（1）））/（1-cos（p
（1）））;
p（3）=p
（2）-（p
（2）-sin（p
（2）））/（1-cos（p
（2）））;
p1=p
（1）-（p
（2）-p
（1））^2/（p（3）-2*p
（2）+p
（1））;
f0=p1-sin（p1）;
ifabs（f0）<
Tol
n=n+1;
p
（1）=p1;
disp（n）近似解以及迭代次数：
假设牛顿迭代法每步的解为pa，重根法每步迭代的解为pb,则
…
21
pa
0.6551
0.4336
0.2881
1.94x10-4
pb
-0.0346
1.38x10-6
3.54x10-11
而采用Steffensen’smethod加速Newton法收敛，可直接一步得到近似解。
因此重根法和Steffensen’smethod可以加速收敛，而Steffensen’smethod的加速效果更明显。
上机作业2
分别用Jacobi、Seidel、Sor（ω=0.8,1.1,1.2,1.3,1.4,1.5）方法求解方程组Ax=b,这里
A=
，b=
x0=[0…0]T;
e=10-5
1、Jacobi法：
迭代格式xi（k+1）=
（bi-
）,i=1,2,…,10,k=1,2,…
>
A=-14*eye（10）+ones（10,10）;
b=-1*ones（10,1）;
x0=zeros（10,1）;
e=10^（-5）;
x=x0;
fori=1:
10
x（i）=（b（i）-A（i,:
）*x0）/A（i,i）+x0（i）;
ifnorm（x-x0）<
e
disp（x）;
disp（k）
return;
x0=x;
ifk>
disp（'
发散'
）
end结果：
2、Seidel法：
）*x）/A（i,i）+x（i）;
end结果：
近似解以及迭代次数
3、Sor迭代法：
迭代格式xi（k+1）=（1-ω）xi（k）+
）,i=1,2,…,10,k=1,2,…
（ω=0.8时）
w=0.8;
x（i）=w*（b（i）-A（i,:
end
结果：
当ω分别取1.1,1.2,1.3,1.4,1.5时，程序不变，只是改变赋值条件中w的值。
当ω=1.1时，ω=1.2时，ω=1.3时，ω=1.4时，ω=1.5时
Jacobi法的迭代列表
k
8
9
x1（k）
0.0769
x2（k）
0.1302
x3（k）
0.167
x4（k）
0.1926
一直计算到29步得到：
x（29）=（0.25,…,0.25）T。
Seidel法的迭代列表如下所示
0.0828
0.0892
0.0961
0.1035
0.1114
0.12
0.1292
0.1392
0.1499
0.1555
0.1611
0.1666
0.172
0.1773
0.1824
0.1872
0.1916
0.1957
0.1992
0.2025
0.2057
0.2087
0.2116
0.2142
0.2166
0.2189
0.221
0.223
0.2248
0.2265
0.2281
0.2296
0.231
0.2323
0.2335
0.2346
0.2356
0.2366
0.2375
一直计算到17步得到：
x（17）=（0.25,…,0.25）T。
Sor迭代法的收敛速度取决于ω的值，取值恰当，便可以加速收敛，在本题中当ω=1.2时，收敛速度最快，其迭代列表如下
0.0923
0.1008
0.1101
0.1203
0.1314
0.1435
0.1568
0.1713
0.1871
0.2043
0.1962
0.2033
0.2101
0.2163
0.2219
0.2267
0.2305
0.2331
0.2342
0.2386
0.2404
0.2419
0.243
0.2438
0.2444
0.2449
0.2455
0.2463
0.2476
0.2475
0.2478
0.248
0.2483
0.2485
0.2488
0.249
0.2492
0.2493
一直计算到9步得到：
x（19）=（0.25,…,0.25）T。
比较以上三种方法可以看出Jacobi的收敛速度慢于Seidel，而Sor则在ω的值选取合适的时候，才能加速收敛。
上机作业3
1、用Newton法与最速下降法求方程组
x2+2y2-2=0
x2=y在（0.8,0.7）附近的根。
解:
1、Newton法：
原理：
由题目可得F（x）=
F′（x）=
应用牛顿迭代法
求解得到近似解。
x0=[0.8;
0.7];
y=-[2*x0
（1）,4*x0
（2）;
2*x0
（1）,-1]\[x0
（1）^2+2*x0
（2）^2-2;
x0
（1）^2-x0
（2）];
x=x0+y;
ifnorm（y）<
=N
2、最速下降法：
原理:
设Ф（x）=f12（x）+f22（x）=（x2+2y2-2）2+（x2-y）2，
（1）求解出其在x（0）点的负梯度方向，G0=grad（x（0））=-（
）T=
-[4x（0）1*（x（0）12+2x（0）22-2）+2x（0）1*（x（0）12-x（0）2）;
8x（0）2*（x（0）12+2x（0）22-2）-2（x（0）12-x（0）2）];
（2）然后求正数λ0,使得Ф（x（0）+λ0G0）=
minФ（x（0）+λ0G0）（λ0>
0）;
（3）取x
（1）=x（0）+λ0G0,则满足Φ（x
（1））<
Φ（x（0））;
（4）x
（1）⇒x（0）,重复上述过程，直到满足Φ（x（k））<
ε。
syma
Q=-[4*x0
（1）*（x0
（1）^2+2*x0
（2）^2-2）+2*x0
（1）*（x0
（1）^2-x0
（2））;
8*x0
（2）*（x0
（1）^2+2*x0
（2）^2-2）-2*（x0
（1）^2-x0
（2））];
y=inline（'
（x1^2+2*x2^2-2）^2+（x1^2-x2）^2'
'
x1'
x2'
）;
u=fminbnd（@（a）y（x0
（1）+a*Q
（1）,x0
（2）+a*Q
（2））,0,100）;
x=x0+u*Q;
if（y（x
（1）,x
（2））<
e）
disp（k）;
上机作业4
1.用幂法与反幂法求矩阵A的按模最大、最小特征值与对应的特征向量。
a）幂法（求模最大特征值与对应特征向量）
A=[4,1,1,1;
1,3,-1,1;
1,-1,2,0;
1,1,0,2];
v=[1;
1;
1];
eps=10^（-5）;
lamda=0;
err=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;
disp（lamda）;
disp（v）;
disp（k）
模最大的特征值和对应的特征向量及迭代次数
b）反幂法（求模最小特征值与对应特征向量）
v0=[1;
[m,i]=max（abs（v0））;
lam0=m;
u0=v0/lam0;
v1=A\u0;
[m1,i]=max（abs（v1））;
lam1=m1;
u1=v1/lam1;
ifabs（1/lam1-1/lam0）<
eps
break
u0=u1;
v0=v1;
lam0=lam1;
lam=1/lam1;
disp（lam）;
disp（u0）;
disp（k）
模最大的特征值和对应的特征向量及迭代次数
2.用Householder变换求矩阵A的QR分解，并用QR方法做3次迭代
Q=eye（4）;
c=norm（A（k:
4,k））;
w=zeros（4,1）;
w（k+1:
4）=A（k+1:
4,k）;
w（k）=A（k,k）+sign（A（k,k））*c;
P=eye（4）-2*w*w'
/norm（w）^2;
Q=Q*P;
A=P*A;
disp（Q）;
disp（A）
矩阵Q以及R（A）
QR法做3次迭代
程序以及运算结果
=3
A=A*Q;
disp（A）
迭代3次后的结果：
上机作业5
目的：
观察lagrange插值的Runge现象。
对于函数f（x）=
进行lagrange插值。
取不同的等分数n（n=5,n=10），将区间[-1，1]n等分，取等距节点。
把插值多项式的曲线画在同一张图上进行比较。
x0=[-1:
0.4:
y0=1./（1+25*x0.^2）;
x=[-1:
0.02:
n=length（x0）;
m=length（x）;
m
z=x（i）;
s=0;
n
t=1;
forj=1:
ifj~=k
t=t*（z-x0（j））/（x0（k）-x0（j））;
s=t*y0（k）+s;
y（i）=s;
plot（x,y,'
-'
holdon
xx0=[-1:
0.2:
yy0=1./（1+25*xx0.^2）;
x1=[-1:
n=length（xx0）;
m=length（x1）;
z=x1（i）;
s=0;
t=t*（z-xx0（j））/（xx0（k）-xx0（j））;
s=t*yy0（k）+s;
y1（i）=s;
plot（x1,y1,'
-.'
运行结果
上机作业6
掌握数值求积公式。
用Romberg求积公式计算积分
，使其精度达到10-4。
a=1;
b=3;
eps=10^（-4）;
M=1;
h=b-a;
J=0;
R=zeros（4,4）;
R（1,1）=h*（exp（a）*cos（a）+exp（b）*cos（b））/2;
while（err>
J=J+1;
h=h/2;
x=a+h:
2*h:
b-h;
R（J+1,1）=R（J,1）/2+h*sum（exp（x））;
forK=1:
min（3,J）
R（J+1,K+1）=R（J+1,K）+（R（J+1,K）-R（J,K））/（4^K-1）;
if（J>
3）
err=abs（R（J+1,4）-R（J,4））;
quad=R（J,4）;
disp（quad）;
disp（R）
积分近似值以及Romberg表：
升级会员 