数值分析上机作业.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)
再用求重根的方法求零点
1.
Newton's method
原理:
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.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+