数值分析上机作业.docx
- 文档编号:25032855
- 上传时间:2023-06-04
- 格式:DOCX
- 页数:60
- 大小:126.51KB
数值分析上机作业.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=x0 a)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+
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 分析 上机 作业