数值分析上机作业.docx
《数值分析上机作业.docx》由会员分享，可在线阅读，更多相关《数值分析上机作业.docx（60页珍藏版）》请在冰点文库上搜索。
数值分析上机作业
数值分析上机实验报告
姓名：
班级：
学号：
第一次上机：
1.题目：
分别用不动点法和牛顿法求解方程
x-exp（x）+4=0
Fixedpointconvergetonegativeroot
Theprincipal：
anumberpisfixedpointforagivenfunctiongifg（p）=p.AndwecanchangeTheequationtofixedpointform.asforfixedpointif
pkmakeg（pk）=pk.
Wecanchangetheformto:
g1=exp（g0）-4
usetheinitialnumberas1
Thecode:
g0=1;
g1=exp（g0）-4;
whileabs（g1-g0）>0.00001
g0=g1;
g1=exp（g0）-4;
end
g1
theresults:
>>w1_fixed_point
g1=
-3.9813
结果分析
通过控制迭代式的形式可以决定收敛到正跟或负根。
Fixedpointconvergetopositiveroot
Theprincipal：
Wecanchangeanewfixedpointformtoproduceapositiveroot
Wecanchangetheformto:
g1=log（4+g0）;
useinitialnumberas-3.9
Thecode:
g0=-3.9;
g1=log（4+g0）;
whileabs（g1-g0）>0.00001
g0=g1;
g1=log（g0+4）;
end
g1
theresults:
>>w1_fixed_point_positive
g1=
1.7490
Newton’smethodtofindnegativeroots
Theprincipal：
WederiveNewton’smethod,basedontaylorpolynomials.newton’smethodisakindoffixedpointmethodwhichhastwoorderconvergence.
Wecanchangetheformto:
g1=g0-（g0-exp（g0）+4）/（1-exp（g0））
useinitialnuberas-3
Thecode:
g0=-3;
g1=g0-（g0-exp（g0）+4）/（1-exp（g0））;
whileabs（g1-g0）>0.00001
g0=g1;
g1=g0-（g0-exp（g0）+4）/（1-exp（g0））
end
theresults:
>>w1_newton
g1=
-3.981342639636226
g1=
-3.981339370911418
结果分析：
收敛速度明显快于不动点（结果太长没打出）。
Newton’smethodtofindpositiveroots
Theprincipal：
Useinitialnumberas1
Thecode:
g0=1;
g1=g0-（g0-exp（g0）+4）/（1-exp（g0））;
whileabs（g1-g0）>0.00001
g0=g1;
g1=g0-（g0-exp（g0）+4）/（1-exp（g0））;
end
g1
theresults:
>>w1_newton
g1=
1.749031386012702
Multipleroots
结果分析：
通过确定初始值得正负来控制收敛到哪个根。
2.p100#9题
题目：
Useeachofthefollowingmethodstofindasolutionin[0.1,1]accuratetowithin10^-4for
600*x^4-550*x^3+200*x^2-20*x-1=0
a.bisectionb.newton’smethod
c.secantmethodd.methodoffalseposition
e.mullers’smethed
a.Bisection
原理：
Themethodcallsforarepeatedhalvingofsubintervalsof[a,b]and,artificialeachstep,locatingthehalfcontainingtherootofequation.
代码：
Thecode:
%bisection
functionanser=Bisection（f,a,b）
tol=1.0e-3;
fa=subs（f,a）;
fmeans=subs（f,（a+b）/2）;
if（fa*fmeans>0）
t=（a+b）/2;
anser=Bisection（f,t,b）;
else
if（fa*fmeans==0）
anser=（a+b）/2;
else
if（abs（b-a）<=tol）
anser=（b+3*a）/4;
else
s=（a+b）/2;
anser=Bisection（f,a,s）;
end
end
end
运行结果：
>>Bisection（'600*x^4-550*x^3+200*x^2-20*x-1',-5,5）
ans=
.0358********
b.newton’smethod
原理:
X=（600*x^4-550*x^3+200*x^2-1）/20
Useinitialnumberas1
代码：
g0=1;
g1=（600*g0^4-550*g0^3+200*g0^2-1）/20;
whileabs（g1-g0）>0.00001
g0=g1;
g1=g0-（g0-exp（g0）+4）/（1-exp（g0））;
end
g1
结果：
>>w1_newton
g1=
1.749031386012702
c.secant
原理:
Useinterval[-5,5]
代码：
functionanser=Secant（f,a,b）
tol=1.0e-5;
error=0.1;
fa=subs（sym（f）,a）;
fb=subs（sym（f）,b）;
anser=a-（b-a）*fa/（fb-fa）;
while（error>tol）
x1=anser;
fx=subs（sym（f）,x1）;
s=fx*fa;
if（s>0）
anser=b-（x1-b）*fb/（fx-fb）;
else
anser=a-（x1-a）*fa/（fx-fa）;
end
error=abs（anser-x1）;
end
end
结果：
Secant（'600*x^4-550*x^3+200*x^2-20*x-1',-5,5）
ans=
0.277662174145539
d.methodoffalseposition
原理:
Useinterval[-5,5]
代码：
%methodofFalsePosition
functionanser=Position（f,a,b）
tol=1.0e-5;
f1=subs（sym（f）,findsym（sym（f））,a）;
f2=subs（sym（f）,findsym（sym（f））,b）;
if（f1==0）
anser=a;
end
if（f2==0）
anser=b;
end
tol1=1;
r1=a;
r2=b;
fv=subs（sym（f）,findsym（sym（f））,a）;
while（tol1>tol）
f2=subs（sym（f）,findsym（sym（f））,r2）;
anser=r2-（r2-r1）*f2/（f2-fv）;
fr=subs（sym（f）,findsym（sym（f））,anser）;
if（f2*fr<0）
tol1=abs（anser-r2）;
r1=r2;
r2=anser;
fv=subs（sym（f）,findsym（sym（f））,r1）;
else
tol1=abs（anser-r2）;
r2=anser;
fv=0.5*subs（sym（f）,findsym（sym（f））,r1）;
end
end
结果：
>>Position（'600*x^4-550*x^3+200*x^2-20*x-1',-5,5）
ans=
0.277662174145539
结果分析：
牛顿法的收敛速度一般较快。
3.
题目：
应用newton法求f（x）的零点，e=10^-6
f（x）=x-sin（x）
再用求重根的方法求零点
１．
Newton＇ｓ ｍｅｔｈｏｄ
原理:
Theiterativeequation
g1=g0-（g0-sin（g0））/（1-cos（g0））;
use1asinitialnumber
代码：
clear;
g0=1;
g1=g0-（g0-sin（g0））/（1-cos（g0））;
whileabs（g1-g0）>0.000001
g0=g1;
g1=g0-（g0-sin（g0））/（1-cos（g0））
end
结果：
>>w3_newton
g1=
1.4990e-006
结果分析：
牛顿法在重根是为一介收敛，失去速度优势。
２．Mutiplerootsmethods1
原理:
g1=g0-（g0-sin（g0））/（1-cos（g0））;
becausegisazerooff（x）ofmultiplicity3
sotheiterationequationis
g1=g0-3*（g0-sin（g0））/（1-cos（g0））
use1asinitialnumber
代码：
clear;
g0=1;
g1=g0-（g0-sin（g0））/（1-cos（g0））;
whileabs（g1-g0）>0.000001
g0=g1;
g1=g0-3*（g0-sin（g0））/（1-cos（g0））
end
结果：
>>w3_mutiple_root
g1=
-0.0095
g1=
2.8751e-008
g1=
6.3996e-009
3．Mutiplerootsmethods2
原理:
g（x）=x–u（x）/u’（x）
=x–f（x）*f’（x）/f’（x）^2–f（x）*f’’（x）
代码：
functionanser=multiple（f,x0）
tol=1.0e-5;
df=diff（sym（f））;
df2=diff（df）;
x1=x0;
error=0.1;
while（error>tol）
fx=subs（f,x1）;
df=subs（df,x1）;
df2=subs（df2,x1）;
anser=x1-fx*df/（df^2-fx*df2）;
error=abs（anser-x1）;
x1=anser;
end
end
结果：
>>Newton_Tay（'x-sin（x）',1）
ans=
0.0302
结果分析：
两种方法的收敛速度都快于牛顿法。
重根公式有效。
4.
题目：
应用newton法求f（x）的零点，e=10^-6这里f（x）=x–sin（x）
再用steffensen’smethod加速其收敛
（1）newton法
原理:
Theiterativeequation
g1=g0-（g0-sin（g0））/（1-cos（g0））;
use1asinitialnumber
代码：
clear;
g0=1;
g1=g0-（g0-sin（g0））/（1-cos（g0））;
whileabs（g1-g0）>0.000001
g0=g1;
g1=g0-（g0-sin（g0））/（1-cos（g0））
end
结果：
>>w3_newton
g1=
1.4990e-006
（2）steffensen’smethod
原理:
Pn={Δ2}（Pn）
代码：
p1=1;
p0=-1;
n=0;
whilen<100&&abs（p1-p0）>0.000001
p1=p0-sin（p0）;
p2=p1-sin（p1）;
p0=p0-（p1-p0）^2/（p2-2*p1+p0）;
n=n+1;
end
p0
结果：
>>w4_steffensen
p0=
0.0358
p0=
-1.6324e-009
p0=
-2.0680e-025
结果分析：
通过加速收敛速度明显加快。
第二次上机
1.
题目：
观察lagrange差值的runge现象
用Neville’s迭代差值算法，对于函数
f（x）=1/（1+25*x^2）,-1<=x<=1
进行lagrange插值。
取不同的等分数
n=5，10，将区间[-1,1]n等份，却等距节点。
把f（x）和插值多项式的曲线画在同一张图上进行比较。
当n=5
原理：
n=5
代码：
xx=-1:
0.01:
1;
yy=zeros（1,201）;
forii=1:
201
n=6;
Q=zeros（n,n）;
x=-1:
（2/（n-1））:
1
x0=xx（ii）;
fori=1:
n
Q（i,1）=1/（1+25*x（i）*x（i））;
end
fori=2:
n
forj=2:
i
Q（i,j）=（（x0-x（i-j+1））*Q（i,j-1）-（x0-x（i））*Q（i-1,j-1））/（x（i）-x（i-j+1））;
end
end
yy（ii）=Q（n,n）;
end
y=1./（1+25.*xx.*xx）;
plot（[xx,xx],[yy,y]）,gridon
结果：
结果分析：
当阶数较少（5阶时）runge现象并不明显。
当n=10
原理：
n=10
代码：
xx=-1:
0.01:
1;
yy=zeros（1,201）;
forii=1:
201
n=11;
Q=zeros（n,n）;
x=-1:
（2/（n-1））:
1
x0=xx（ii）;
fori=1:
n
Q（i,1）=1/（1+25*x（i）*x（i））;
end
fori=2:
n
forj=2:
i
Q（i,j）=（（x0-x（i-j+1））*Q（i,j-1）-（x0-x（i））*Q（i-1,j-1））/（x（i）-x（i-j+1））;
end
end
yy（ii）=Q（n,n）;
end
y=1./（1+25.*xx.*xx）;
plot（[xx,xx],[yy,y]）,gridon
结果：
结果分析：
阶数为10时观察到明显的runge现象，在边缘部分误差突然增大。
2.
题目p155#28题
Theupperportionofthisnoblebeastistobeapproximatedusingclampedcubicsplineinterpolants.thecurveisdrawnonagridfromwhichthetableisconstructed.usealgorithm3.5toconstructtheclampedcubicsplines.
原理：
Givenafunctionfdefinedon[a,b]andasetofnodesa=x0a）S（x）isacubicpolynomial,denotedSi（x）,onthe
subinterval[xi,xi+1]foreachi=0,1,…,n-1;
代码：
function[a,b,c,d]=Clampcubicspline（x,y,FPO,FPN）
n=length（x）;
m=length（y）;
ifm~=n
end
n=n-1;
a=y;
fori=1:
n
h（i）=x（i+1）-x（i）;%用h（i）记录x的差值
end
Alph=zeros（1,n+1）;
Alph
（1）=3*（a
（2）-a
（1））/h
（1）-3*FPO;
Alph（n+1）=3*FPN-3*（a（n+1）-a（n））/h（n）;
fori=2:
n
Alph（i）=3*（a（i+1）-a（i））/h（i）-3*（a（i）-a（i-1））/h（i-1）;
end
l=zeros（1,n+1）;
u=zeros（1,n）;
z=zeros（1,n+1）;
l
（1）=2*h
（1）;
u
（1）=0.5;
z
（1）=Alph
（1）/l
（1）;
fori=2:
n
l（i）=2*（x（i+1）-x（i-1））-h（i-1）*u（i-1）;
u（i）=h（i）/l（i）;
z（i）=（Alph（i）-h（i-1）*z（i-1））/l（i）;
end
l（n+1）=h（n）*（2-u（n））;
z（n+1）=（Alph（n+1）-h（n）*z（n））/l（n+1）;
c=zeros（1,n+1）;
c（n+1）=z（n+1）;
forj=n:
-1:
1
c（j）=z（j）-u（j）*c（j+1）;
b（j）=（a（j+1）-a（j））/h（j）-h（j）*（c（j+1）+2*c（j））/3;
d（j）=（c（j+1）-c（j））/（3*h（j））;
end
a=a（1:
n）;
c=c（1:
n）;
%curve1
x1=[125678101317];
y1=[3.03.73.94.25.76.67.16.74.5];
FPO=1.0;
FPN=-2/3;
[a,b,c,d]=Clampcubicspline（x1,y1,FPO,FPN）;
n=length（x1）;
fori=1:
n-1
t=（x1（i+1）-x1（i））/50;
xx=x1（i）:
t:
x1（i+1）;
yy=a（i）+b（i）*（xx-x1（i））+c（i）*（xx-x1（i））.^2+d（i）*（xx-x1（i））.^3;
plot（xx,yy,'r'）
holdon
end
%curve2
x2=[17202324252727.7];
y2=[4.57.06.15.65.85.24.1];
FPO=3;
FPN=-4;
[a,b,c,d]=Clampcubicspline（x2,y2,FPO,FPN）;
n=length（x2）;
fori=1:
n-1
t=（x2（i+1）-x2（i））/50;
xx=x2（i）:
t:
x2（i+1）;
yy=a（i）+b（i）*（xx-x2（i））+c（i）*（xx-x2（i））.^2+d（i）*（xx-x2（i））.^3;
plot（xx,yy,'c'）
holdon
end
%curve3
x3=[27.7282930];
y3=[4.14.34.13.0];
FPO=1/3;
FPN=-3/2;
[a,b,c,d]=Clampcubicspline（x3,y3,FPO,FPN）;
n=length（x3）;
fori=1:
n-1
t=（x3（i+1）-x3（i））/50;
xx=x3（i）:
t:
x3（i+1）;
yy=a（i）+b（i）*（xx-x3（i））+c（i）*（xx-x3（i））.^2+d（i）*（xx-x3（i））.^3;
plot（xx,yy,'g'）
holdon
end
plot（x1,y1,'*'）
holdon
plot（x2,y2,'k*'）
holdon
plot（x3,y3,'m*'）
title（'clampedcubicsplineinterpolants近似狗的背部曲线'）
gridon
结果：
结果分析：
结果与原图较为接近。
3.
题目p15510
UseRombergintrgrationtocomputethefollowingapproximationsto
原理：
RombergintegrationusestheCompositeTrapezoidal
ruletogivepreliminaryapproximationsandthen
appliestheRichardsonextrapolationprocessto
improvetheapproximations.
代码：
fun='sqrt（1+（cos（x））.^2）';
a=0,b=48,tol=1.e-8;
k=0;
n=1;
h=b-a;
T=h/2*（sqrt（1+（cos（a））.^2）+sqrt（1+（cos（b））.^2））;
error=1;
whileerror>=tol
k=k+1;
h=h/2;
tmp=0;
fori=1:
n
tmp=tmp+sqrt（1+（cos（a+（2*i-1）*h））.^2）;
end
T（k+1,1）=T（k）/2+h*tmp;
forj=1:
k
T（k+1,j+1）=T（k+1,j）+（T（k+1,j）-T（k,j））/（4^j-1）;
end
n=n*2;
error=abs（T（k+
升级会员 