数值分析计算实习题.docx
- 文档编号:10115513
- 上传时间:2023-05-23
- 格式:DOCX
- 页数:25
- 大小:89.61KB
数值分析计算实习题.docx
《数值分析计算实习题.docx》由会员分享,可在线阅读,更多相关《数值分析计算实习题.docx(25页珍藏版)》请在冰点文库上搜索。
数值分析计算实习题
插值法
1.下列数据点的插值
x01491625364964
y012345678
可以得到平方根函数的近似,在区间[0,64]上作图.
(1)用这9个点作8次多项式插值Ls(x).
(2)用三次样条(第一边界条件)程序求S(x).
从得到结果看在[0,64]上,哪个插值更精确;在区间[0,1]上,两种插值哪个更精确?
解:
(1)拉格朗日插值多项式,求解程序如下
symsxl;
x1=[01491625364964];
y1=[012345678];
n=length(x1);
Ls=sym(0);
fori=1:
n
l=sym(y1(i));
fork=1:
i-1
l=l*(x-x1(k))/(x1(i)-x1(k));
end
fork=i+1:
n
l=l*(x-x1(k))/(x1(i)-x1(k));
end
Ls=Ls+l;
end
Ls=simplify(Ls)%为所求插值多项式Ls(x).
输出结果为
Ls=
-24221063/63504000*x^2+95549/72072*x-1/3048192000*x^8-2168879/435456000*x^4+19/283046400*x^7+657859/10886400*x^3+33983/152409600*x^5-13003/2395008000*x^6
(2)三次样条插值,程序如下
x1=[01491625364964];
y1=[012345678];
x2=[0:
1:
64];
y3=spline(x1,y1,x2);
p=polyfit(x2,y3,3);%得到三次样条拟合函数
S=p
(1)+p
(2)*x+p(3)*x^2+p(4)*x^3%得到S(x)
输出结果为:
S=
2288075067923491/73786976294838206464-2399112304472833/576460752303423488*x+4552380473376713/18014398509481984*x^2+999337332656867/1125899906842624*x^3
(3)在区间[0,64]上,分别对这两种插值和标准函数作图,
plot(x2,sqrt(x2),'b',x2,y2,'r',x2,y3,'y')
蓝色曲线为y=
函数曲线,红色曲线为拉格朗日插值函数曲线,黄色曲线为三次样条插值曲线
可以看到蓝色曲线与黄色曲线几乎重合,因此在区间[0,64]上三次样条插值更精确。
在[0,1]区间上由上图看不出差别,不妨代入几组数据进行比较,取x4=[0:
0.2:
1]
x4=[0:
0.2:
1];
sqrt(x4)%准确值
subs(Ls,'x',x4)%拉格朗日插值
spline(x1,y1,x4)%三次样条插值
运行结果为
ans=
00.44720.63250.77460.89441.0000
ans=
00.25040.47300.67060.84551.0000
ans=
00.24290.46300.66170.84031.0000
从这几组数值上可以看出在[0,1]区间上,拉格朗日插值更精确。
数据拟合和最佳平方逼近
2.有实验给出数据表
x0.00.10.20.30.50.81.0
y1.00.410.500.610.912.022.46
试求3次、4次多项式的曲线拟合,再根据数据曲线形状,求一个另外函数的拟合曲线,用图示数据曲线及相应的三种拟合曲线。
解:
(1)三次拟合,程序如下
symx;
x0=[0.00.10.20.30.50.81.0];
y0=[1.00.410.500.610.912.022.46];
cc=polyfit(x0,y0,3);
S3=cc
(1)+cc
(2)*x+cc(3)*x^2+cc(4)*x^3%三次拟合多项式
xx=x0
(1):
0.1:
x0(length(x0));
yy=polyval(cc,xx);
plot(xx,yy,'--');
holdon;
plot(x0,y0,'x');
xlabel('x');
ylabel('y');
运行结果
S3=
-7455778416425083/1125899906842624+1803512222945437/140737488355328*x-655705280524945/140737488355328*x^2+4172976892199509/4503599627370496*x^3
图像如下
(2)4次多项式拟合
symx;
x0=[0.00.10.20.30.50.81.0];
y0=[1.00.410.500.610.912.022.46];
cc=polyfit(x0,y0,4);
S3=cc
(1)+cc
(2)*x+cc(3)*x^2+cc(4)*x^3+cc(5)*x^4
xx=x0
(1):
0.1:
x0(length(x0));
yy=polyval(cc,xx);
plot(xx,yy,'r');
holdon;
plot(x0,y0,'x');
xlabel('x');
ylabel('y');
运行结果
S3=
3248542900396215/1125899906842624-3471944732519151/281474976710656*x+4580931990070637/281474976710656*x^2-5965836931688425/1125899906842624*x^3+8491275464650307/9007199254740992*x^4
图像如下
(3)另一个拟合曲线,
新建一个M-file
程序如下:
function[C,L]=lagran(x,y)
w=length(x);
n=w-1;
L=zeros(w,w);
fork=1:
n+1
V=1;
forj=1:
n+1
ifk~=j
V=conv(V,poly(x(j)))/(x(k)-x(j));
end
end
L(k,:
)=V;
end
C=y*L
在命令窗口中输入以下的命令:
x=[0.00.10.20.30.50.81.0];
y=[1.00.410.500.610.912.022.46];
cc=polyfit(x,y,4);
xx=x
(1):
0.1:
x(length(x));
yy=polyval(cc,xx);
plot(xx,yy,'r');
holdon;
plot(x,y,'x');
xlabel('x');
ylabel('y');
x=[0.00.10.20.30.50.81.0];
y=[0.940.580.470.521.002.002.46];%y中的值是根据上面两种拟合曲线而得到的估计数据,不是真实数据
[C,L]=lagran(x,y);
xx=0:
0.01:
1.0;
yy=polyval(C,xx);
holdon;
plot(xx,yy,'b',x,y,'.');
图像如下
解线性方程组的直接解法
3.线性方程组Ax=b的A及b为
A=
,b=
,则解x=
.用MATLAB内部函数求detA及A的所有特征值和cond(A)2.若令
A+δA=
,求解(A+δA)(x+δx)=b,输出向量x和||δx||2.从理论结果和实际计算两方面分析线性方程组Ax=b解得相对误差||δx||2/||x||2及A的相对误差||δA||2/||A||2的关系.
解:
(1)程序如下
clear;
A=[10787;7565;86109;75910];
det(A)
cond(A,2)
eig(A)
输出结果为
ans=
1
ans=
2.9841e+003
ans=
0.0102
0.8431
3.8581
30.2887
(2)程序如下
A=[1078.17.2;7.085.0465;85.989.899;6.99599.98];
b=[32233331]';
x0=[1111]';
x=A\b%扰动后方程组的解
x1=x-x0%δx的值
norm(x1,2)%δx的2-范数
运行结果为
x=
-9.5863
18.3741
-3.2258
3.5240
x1=
-10.5863
17.3741
-4.2258
2.5240
ans=
20.9322
程序如下
A=[1078.17.2;7.085.0465;85.989.899;6.99599.98];
A0=[10787;7565;86109;75910];
b=[32233331]';
x0=[1111]';
x=A\b;
x1=x-x0;
norm(x1,2);
A1=A-A0;%δA的值
norm(x1,2)/norm(x0,2)%||δx||2/||x||2的值
norm(A1,2)/norm(A0,2)%||δA||2/||A||2的值
输出结果为
ans=
10.4661
ans=
0.0076
可见A相对误差只为0.0076,而得到的结果x的相对误差就达到了10.4661,该方程组是病态的,A的条件数为2984.1远远大于1,当A只有很小的误差就会给结果带来很大的影响。
非线性方程数值解法
4.求下列方程的实根
(1)x^2-3x+2-e^x=0;
(2)x^3+2x^2+10x-20=0.
要求:
(1)设计一种不动点迭代法,要使迭代序列收敛,然后再用斯特芬森加速迭代,计算到|x(k)-x(k-1)|<10^(-8)为止。
(2)用牛顿迭代,同样计算到|x(k)-x(k-1)|<10^(-8)。
输出迭代初值及各次迭代值和迭代次数k,比较方法的优劣。
解:
(1)先用画图的方法估计根的范围
ezplot('x^2-3*x+2-exp(x)');
gridon;
可以估计到方程的根在区间(0,1);选取迭代初值为x0=0.5;
构造不动点迭代公式x(k+1)=(x(k)^2+2-e^x(k))/3;
ψ(x)=(x^2+2-e^x)/3;
当0 程序如下: formatlong; f=inline('(x^2+2-exp(x))/3') disp('x='); x=feval(f,0.5); disp(x); Eps=1E-8; i=1; while1 x0=x; i=i+1; x=feval(f,x); disp(x); ifabs(x-x0) break; end end i,x 运行结果为 f= Inlinefunction: f(x)=(x^2+2-exp(x))/3 x= 0.20042624309996 0.27274906509837 0.25360715658413 0.25855037626494 0.25726563633509 0.25759898516219 0.25751245451483 0.25753491361525 0.25752908416796 0.25753059723833 0.25753020451046 0.25753030644564 0.25753027998767 0.25753028685501 i= 14 x= 0.25753028685501 斯特芬森加速法,程序如下: formatlong; f=inline('x-((x^2+2-exp(x))/3-x)^2/((((x^2+2-exp(x))/3)^2+2-exp((x^2+2-exp(x))/3))/3-2*(x^2+2-exp(x))/3+x)'); disp('x='); x=feval(f,0.5); disp(x); Eps=1E-8; i=1; while1 x0=x; i=i+1; x=feval(f,x); disp(x); ifabs(x-x0) break; end end i,x 运行结果为x= 0.25868442756579 0.25753031771981 0.25753028543986 0.25753028543986 i= 4 x= 0.25753028543986 牛顿迭代法,程序如下: formatlong; x=sym('x'); f=sym('x^2-3*x+2-exp(x)'); df=diff(f,x); FX=x-f/df; Fx=inline(FX); disp('x='); x1=0.5; disp(x1); Eps=1E-8; i=0; while1 x0=x1; i=i+1; x1=feval(Fx,x1); disp(x1); ifabs(x1-x0) break; end end i,x1 运行结果如下: x= 0.50000000000000 0.25368870241829 0.25752890079471 0.25753028543968 0.25753028543986 i= 4 x1= 0.25753028543986 (2)先用画图的方法估计根的范围 ezplot('x^3+2*x^2+10*x-20'); gridon; 根大约在区间(1,2);选取初值x0=1.5; 构造不动点迭代公式x(k+1)=(-2x(k)^2-10x(k)+20)^1/3; ψ(x)=(-2x^2-10x+20)^1/3; 程序如下: formatlong; f=inline('(-2*x^2-10*x+20)^1/3') disp('x='); x=feval(f,1.5); disp(x) Eps=1E-8; i=1; while1 x0=x; i=i+1; x=feval(f,x); disp(x); ifabs(x-x0)>1E10 break; end ifabs(x-x0) break; end end i,x 运行结果: f= Inlinefunction: f(x)=(-2*x^2-10*x+20)^1/3 x= 0.166********667 6.09259259259259 -38.38843164151806 -8.478196837919431e+002 -4.763660785374071e+005 -1.512815059604763e+011 i= 6 x= -1.512815059604763e+011 迭代6次后x的值大得令人吃惊,表明构造的式子并不收敛. 也无法构造出收敛的不动点公式 牛顿迭代法,程序如下: formatlong; x=sym('x'); f=sym('x^3+2*x^2+10*x-20'); df=diff(f,x); FX=x-f/df; Fx=inline(FX); disp('x='); x1=0.5; disp(x1); Eps=1E-8; i=0; while1 x0=x1; i=i+1; x1=feval(Fx,x1); disp(x1); ifabs(x1-x0) break; end end i,x1 运行结果: x= 1.50000000000000 1.37362637362637 1.36881481962396 1.36880810783441 1.36880810782137 i= 4 x1= 1.36880810782137 比较三种方法,牛顿法的收敛性比较好,相比不动点迭代法要构造出收敛的公式比较难,牛顿法迭代次数也较少,收敛速度快,只是对初值的要求很高,几种方法各有利弊,具体采用哪种也需因题而异。 常微分方程初值问题数值解法 5.给定初值问题 y’=-50y+50x^2+2x, ; y(0)=1/3; 用经典的四阶R-K方法解该问题,步长分别取h=0.1,0.025,0.01,计算并打印x=0.1i(i=0,1,…,10)各点的值,与准确值y(x)=1/3e^(-50x)+x^2比较。 解: 取步长h=0.1,程序如下: %经典的四阶R-K方法 clear; F='-50*y+50*x^2+2*x'; a=0;b=1; h=0.1; n=(b-a)/0.1; X=a: 0.1: b; Y=zeros(1,n+1); Y (1)=1/3; fori=1: n x=X(i); y=Y(i); K1=h*eval(F); x=x+h/2; y=y+K1/2; K2=h*eval(F); y=Y(i)+K2/2; K3=h*eval(F); x=X(i)+h; y=Y(i)+K3; K4=h*eval(F); Y(i+1)=Y(i)+(K1+2*K2+2*K3+K4)/6; end %准确值 temp=[]; f=dsolve('Dy=-50*y+50*x^2+2*x','y(0)=1/3','x'); df=zeros(1,n+1); fori=1: n+1 temp=subs(f,'x',X(i)); df(i)=double(vpa(temp)); end disp('步长四阶经典R-K法准确值'); disp([X',Y',df']); 运行结果: 步长四阶经典R-K法准确值 1.0e+010* 00.000000000033330.00000000003333 0.000000000010000.000000000460550.00000000000122 0.000000000020000.000000006306250.00000000000400 0.000000000030000.000000086404940.00000000000900 0.000000000040000.000001184363000.00000000001600 0.000000000050000.000016235451100.00000000002500 0.000000000060000.000222560671340.00000000003600 0.000000000070000.003050935427780.00000000004900 0.000000000080000.041823239217400.00000000006400 0.000000000090000.573326903478090.00000000008100 0.000000000100007.859356300837710.00000000010000 %画图观察结果 figure; plot(X,df,'k*',X,Y,'--r'); gridon; title('四阶经典R-K法解常微分方程'); legend('准确值','四阶经典R-K法'); 当x值接近1的时候,偏离准确值太大。 当步长h=0.025时,将上面程序中的h改为0.025即可,运行结果: 步长四阶经典R-K法准确值 00.333333333333330.33333333333333 0.100000000000000.103135240342880.01224598233303 0.200000000000000.044285273275990.04001513330992 0.300000000000000.051967957555070.09000010196744 0.400000000000000.093957311494390.16000000068705 0.500000000000000.160345314357570.25000000000463 0.600000000000000.248085701305570.36000000000003 0.700000000000000.356241884727580.49000000000000 0.800000000000000.484525906616270.64000000000000 0.900000000000000.632849233007510.81000000000000 1.000000000000000.801184643742061.00000000000000 图像如下: 当步长h=0.01时,将上面程序中的h改为0.01即可,运行结果: 步长四阶经典R-K法准确值 00.333333333333330.33333333333333 0.100000000000000.202357204861110.01224598233303 0.200000000000000.128817001907910.04001513330992 0.300000000000000.097991826678500.09000010196744 0.400000000000000.100949467750240.16000000068705 0.500000000000000.132********4700.25000000000463 0.600000000000000.188********1990.36000000000003 0.700000000000000.268139302784310.49000000000000 0.800000000000000.369481660283190.64000000000000 0.900000000000000.491957621994750.81000000000000 1.000000000000000.635121421679101.00000000000000
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 分析 计算 实习