数值分析计算实习题.docx
《数值分析计算实习题.docx》由会员分享,可在线阅读,更多相关《数值分析计算实习题.docx(28页珍藏版)》请在冰豆网上搜索。
数值分析计算实习题
插值法
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
七年级英语期末考试质量分析
一、试卷分析:
本次试卷的难易程度定位在面向大多数学生。
该份试卷紧扣教材,突出重点,注重对基础知识