数值分析计算实习题.docx
- 文档编号:9273188
- 上传时间:2023-02-03
- 格式:DOCX
- 页数:25
- 大小:86.05KB
数值分析计算实习题.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=
23491/304472833/8*x+76713/*x^2+6867/42624*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=
-25083/42624+45437/5328*x-4945/5328*x^2+99509/70496*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=
96215//0656*x+70637/0656*x^2-88425/42624*x^3+50307/40992*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数
运行结果为
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)先用画图的方法估计根的X围
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.996 0.837 0.413 0.494 0.509 0.219 0.483 0.525 0.796 0.833 0.046 0.564 0.767 0.501 i= 14 x= 0.501 斯特芬森加速法,程序如下: 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.579 0.981 0.986 0.986 i= 4 x= 0.986 牛顿迭代法,程序如下: 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.000 0.829 0.471 0.968 0.986 i= 4 x1= 0.986 (2)先用画图的方法估计根的X围 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.667 6.259 -38.806 -8.9431e+002 -4.4071e+005 -1.4763e+011 i= 6 x= -1.4763e+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.000 1.637 1.396 1.441 1.137 i= 4 x1= 1.137 比较三种方法,牛顿法的收敛性比较好,相比不动点迭代法要构造出收敛的公式比较难,牛顿法迭代次数也较少,收敛速度快,只是对初值的要求很高,几种方法各有利弊,具体采用哪种也需因题而异。 常微分方程初值问题数值解法 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.3330.333 0.0000.0550.122 0.0000.6250.400 0.0000.4940.900 0.0000.3000.600 0.0000.1100.500 0.0000.1340.600 0.0000.7780.900 0.0000.7400.400 0.0000.8090.100 0.0007.7710.000 %画图观察结果 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.3330.333 0.0000.2880.303 0.0000.5990.992 0.0000.5070.744 0.0000.4390.705 0.0000.7570.463 0.0000.5570.003 0.0000.7580.000 0.0000.6270.000 0.0000.7510.000 1.0000.2061.000 图像如下: 当步长h=0.01时,将上面程序中的h改为0.01即可,运行结果: 步长四阶经典R-K法准确值 00.3330.333 0.0000.1110.303 0.0000.7910.992 0.0000.8500.744 0.0000.0240.705 0.0000.4700.463 0.0000.1990.003 0.0000.4310.000 0.0000.3190.000 0.0000.4750.000 1.0000.9101.000
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值 分析 计算 实习