数值分析实验报告.docx
《数值分析实验报告.docx》由会员分享,可在线阅读,更多相关《数值分析实验报告.docx(20页珍藏版)》请在冰豆网上搜索。
![数值分析实验报告.docx](https://file1.bdocx.com/fileroot1/2023-6/7/16de49bc-c534-479f-8d3a-f3de7cb87f45/16de49bc-c534-479f-8d3a-f3de7cb87f451.gif)
数值分析实验报告
实验2.1多项式插值的振荡现象
实验目的:
在一个固定的区间上用插值逼近一个函数,显然Lagrange插值中使用的节点越多,插值多项式的次数就越高。
我们自然关心插值多项式的次数增加时,Ln(x)是否也更加靠近被逼近的函数。
Runge给出的一个例子是极著名并富有启发性的。
实验容:
设区间[-1,1]上函数f(x)=1/(1+25x2)。
考虑区间[-1,1]的一个等距划分,分点为xi=-1+2i/n,i=0,1,2,…,n,
则拉格朗日插值多项式为
.
其中,li(x),i=0,1,2,…,n是n次Lagrange插值基函数。
实验步骤与结果分析:
实验源程序
functionChap2Interpolation
%数值实验二:
“实验2.1:
多项式插值的震荡现象”
%输入:
函数式选择,插值结点数
%输出:
拟合函数及原函数的图形
promps={'请选择实验函数,若选f(x),请输入f,若选h(x),请输入h,若选g(x),请输入g:
'};
titles='charpt_2';
result=inputdlg(promps,'charpt2',1,{'f'});
Nb_f=char(result);
if(Nb_f~='f'&Nb_f~='h'&Nb_f~='g')errordlg('实验函数选择错误!
');return;end
result=inputdlg({'请输入插值结点数N:
'},'charpt_2',1,{'10'});
Nd=str2num(char(result));
if(Nd<1)errordlg('结点输入错误!
');return;end
switchNb_f
case'f'
f=inline('1./(1+25*x.^2)');a=-1;b=1;
case'h'
f=inline('x./(1+x.^4)');a=-5;b=5;
case'g'
f=inline('atan(x)');a=-5;b=5;
end
x0=linspace(a,b,Nd+1);y0=feval(f,x0);
x=a:
0.1:
b;y=Lagrange(x0,y0,x);
fplot(f,[ab],'co');
holdon;
plot(x,y,'b--');
xlabel('x');ylabel('y=f(x)oandy=Ln(x)--');
%--------------------------------------------------------------------
functiony=Lagrange(x0,y0,x);
n=length(x0);m=length(x);
fori=1:
m
z=x(i);
s=0.0;
fork=1:
n
p=1.0;
forj=1:
n
if(j~=k)
p=p*(z-x0(j))/(x0(k)-x0(j));
end
end
s=s+p*y0(k);
end
y(i)=s;
end
实验结果分析
(1)增大分点n=2,3,…时,拉格朗日插值函数曲线如图所示。
n=6
n=7
n=8
n=9
n=10
从图中可以看出,随着n的增大,拉格朗日插值函数在x=0附近较好地逼近了原来的函数f(x),但是却在两端x=-1和x=1处出现了很大的振荡现象。
并且,仔细分析图形,可以看出,当n为奇数时,虽然有振荡,但振荡的幅度不算太大,n为偶数时,其振荡幅度变得很大。
通过思考分析,我认为,可能的原因是f(x)本身是偶函数,如果n为奇数,那么Lagrange插值函数Ln(x)的最高次项xn-1是偶次幂,比较符合f(x)本身是偶函数的性质;如果n为偶数,那么Lagrange插值函数Ln(x)的最高次项xn-1是奇次幂,与f(x)本身是偶函数的性质相反,因此振荡可能更剧烈。
(2)将原来的f(x)换为其他函数如h(x)、g(x),结果如图所示。
其中h(x),g(x)均定义在[-5,5]区间上,h(x)=x/(1+x4),g(x)=arctanx。
h(x),n=7
h(x),n=8
h(x),n=9
h(x),n=10
g(x),n=7
g(x),n=8
g(x),n=9
g(x),n=10
分析两个函数的插值图形,可以看出:
随着n的增大,拉格朗日插值函数在x=0附近较好地逼近了原来的函数f(x),但是却在两端x=-5和x=5处出现了很大的振荡现象。
并且,仔细分析图形,可以看出,当n为偶数时,虽然有振荡,但振荡的幅度不算太大,n为奇数时,其振荡幅度变得很大。
原因和上面f(x)的插值类似,h(x)、g(x)本身是奇函数,如果n为偶数,那么Lagrange插值函数Ln(x)的最高次项xn-1是奇次幂,比较符合h(x)、g(x)本身是奇函数的性质;如果n为奇数,那么Lagrange插值函数Ln(x)的最高次项xn-1是偶次幂,与h(x)、g(x)本身是奇函数的性质相反,因此振荡可能更剧烈。
实验3.1多项式最小二乘拟合
实验目的:
编制以函数{xk}k=0,…,n;为基的多项式最小二乘拟合程序。
实验容:
对表中的数据作三次多项式最小二乘拟合。
xi
-1.0
-0.5
0.0
0.5
1.0
1.5
2.0
yi
-4.447
-0.452
0.551
0.048
-0.447
0.549
4.552
取权函数wi≡1,求拟合曲线
中的参数{αk}、平方误差δ2,并作离散据
{xi,yi}的拟合函数的图形。
实验源程序
functionChap3CurveFitting
%数值实验三:
“实验3.1”
%输出:
原函数及求得的相应插值多项式的函数的图像以及参数alph和误差r
x0=-1:
0.5:
2;
y0=[-4.447-0.4520.5510.048-0.4470.5494.552];
n=3;%n为拟合阶次
alph=polyfit(x0,y0,n);
y=polyval(alph,x0);
r=(y0-y)*(y0-y)';%平方误差
x=-1:
0.01:
2;
y=polyval(alph,x);
plot(x,y,'k--');
xlabel('x');ylabel('y0*andpolyfit.y--');
holdon
plot(x0,y0,'*')
gridon;
disp(['平方误差:
',num2str(r)])
disp(['参数alph:
',num2str(alph)])
实验结果
平方误差:
2.1762e-005
参数alph:
1.9991-2.9977-3.9683e-0050.54912
实验4.1
实验目的:
复化求积公式计算定积分.
实验题目:
数值计算下列各式右端定积分的近似值.
实验要求:
(1)若用复化梯形公式、复化Simpson公式和复化Gauss-LegendreI型公式做计算,要求绝对误差限为
,分别利用它们的余项对每种算法做出步长的事前估计.
(2)分别用复化梯形公式,复化Simpson公式和复化Gauss-LegendreI型公式作计算.
(3)将计算结果与精确解做比较,并比较各种算法的计算量.
实验程序:
1.事前估计的Matlab程序如下:
(1).用复化梯形公式进行事前估计的Matlab程序
formatlongg
x=2:
0.01:
3;
f=-4*(3*x.^2+1)./(x.^2-1).^3;%二阶导函数
%plot(x,f)%画出二阶导函数图像
x=2.0;%计算导函数最大值
f=-4*(3*x^2+1)/(x^2-1)^3;
h2=0.5*10^(-7)*12/f;
h=sqrt(abs(h2))%步长
n=1/h;
n=ceil(1/h)+1%选取的点数
formatlongg
x=0:
0.01:
1;
f=8.*(3*x.^2-1)./(x.^2+1).^3;%二阶导函数
%plot(x,f)%画出二阶导函数图像
x=1;%计算导函数最大值
f=8.*(3*x.^2-1)./(x.^2+1).^3;
h2=0.5*10^(-7)*12/f;
h=sqrt(abs(h2))%步长
n=1/h
n=ceil(1/h)+1%选取的点数
formatlongg
x=0:
0.01:
1;
f=log(3).*log(3).*3.^x;%二阶导函数
%plot(x,f);%画出二阶导函数图像
x=1;%计算导函数最大值
f=log(3)*log(3)*3^x;
h2=0.5*10^(-7)*12/f;
h=sqrt(abs(h2))%步长
n=1/h
n=ceil(1/h)+1%选取的点数
formatlongg
x=1:
0.01:
2;
f=2.*exp(x)+x.*exp(x);%二阶导函数
%plot(x,f)%画出二阶导函数图像
x=2;%计算导函数最大值
f=2.*exp(x)+x.*exp(x);
h2=0.5*10^(-7)*12/f;
h=sqrt(abs(h2))%步长
n=1/h
n=ceil(1/h)+1%选取的点数
估计结果步长h及结点数n分别为
h=0.8
n=1793
h=0.6
n=1827
h=0.9
n=2458
h=0.9
n=7020
(2).用复化simpson公式进行事前估计的Matlab程序
formatlongg
x=2:
0.01:
3;
f=-2*((-72*x.^2-24).*(x.^2-1)-192*x.^2.*(x.^2+1))./(x.^2-1).^5;%四阶导函数
x=2.0;
f=-2*((-72*x^2-24)*(x^2-1)-192*x^2*(x^2+1))/(x^2-1)^5;%计算导函数最大值
h4=0.5*10^(-7)*180*16/f;
h=sqrt(sqrt(abs(h4)))%步长
n=1/h;%求分段区间个数
n=2*ceil(1/h)+1%选取的点数
formatlongg
x=0:
0.01:
1;
f=4*((-72*x.^2+24).*(x.^2+1)-192*x.^2.*(-x.^2+1))./(x.^2+1).^5;%四阶导函数
x=1;
f=4*((-72*x^2+24)*(x^2+1)-192*x^2*(-x^2+1))/(x^2+1)^5;%计算导函数最大值
h4=0.5*10^(-7)*180*16/f;
h=sqrt(sqrt(abs(h4)))%步长
n=1/h;%求分段区间个数
n=2*ceil(1/h)+1%选取的点数
formatlongg
x=0:
0.01:
1;
f=log(3)^4*3.^x;%四阶导函数
x=1;
f=log(3)^4*3.^x;%计算导函数最大值
h4=0.5*10^(-7)*180*16/f;
h=sqrt(sqrt(abs(h4)))%步长
n=1/h;%求分段区间个数
n=2*ceil(1/h)+1%选取的点数
formatlongg
x=1:
0.01:
2;
f=4*exp(x)+x.*exp(x);%四阶导函数
plot(x,f)%画出原函数
x=2;
f=4*exp(x)+x.*exp(x);%计算导函数最大值
h4=0.5*10^(-7)*180*16/f;
h=sqrt(sqrt(abs(h4)))
n=1/h;%求分段区间个数
n=2*ceil(1/h)+1%选取的点数
估计结果步长h及结点数n分别为
h=0.13411
n=47
h=0.76542
n=35
h=0.18433
n=29
h=0.18546
n=49
2.积分计算的Matlab程序:
formatlongg
promps={'请选择积分公式,若用复化梯形,请输入T,用复化simpson,输入S,
用复化Gauss_Legendre,输入GL:
'};
result=inputdlg(promps,'charpt4',1,{'T'});
Nb=char(result);
if(Nb~='T'&Nb~='S'&Nb~='GL')
errordlg('积分公式选择错误');
return;
end
result=inputdlg({'请输入积分式题号1-4:
'},'实验4.1',1,{'1'});
Nb_f=str2num(char(result));
if(Nb_f<1|Nb_f>4)
errordlg('没有该积分式');
return;
end
switchNb_f
case1
fun=inline('-2./(x.^2-1)');a=2;b=3;
case2
fun=inline('4./(x.^2+1)');a=0;b=1;
case3
fun=inline('3.^x');a=0;b=1;
case4
fun=inline('x.*exp(x)');a=1;b=2;
end
if(Nb=='T')%用复化梯形公式
promps={'请输入用复化梯形公式应取的步长:
'};
result=inputdlg(promps,'实验4.2',1,{'0.01'});
h=str2num(char(result));
if(h<=0)
errordlg('请输入正确的步长!
');
return;
end
tic;
N=floor((b-a)/h);
detsum=0;
fori=1:
N-1
xk=a+i*h;
detsum=detsum+fun(xk);
end
t=h*(fun(a)+fun(b)+2*detsum)/2;
time=toc;
t
end
if(Nb=='S')%用复化Simpson公式
promps={'请输入用复化Simpson公式应取的步长:
'};
result=inputdlg(promps,'实验4.2',1,{'0.01'});
h=str2num(char(result));
if(h<=0)
errordlg('请输入正确的步长!
');
return;
end
tic;
N=floor((b-a)/h);
detsum_1=0;
detsum_2=0;
fori=1:
N-1
xk_1=a+i*h;
detsum_1=detsum_1+fun(xk_1);
end
fori=1:
N
xk_2=a+h*(2*i-1)/2;
detsum_2=detsum_2+fun(xk_2);
end
t=h*(fun(a)+fun(b)+2*detsum_1+4*detsum_2)/6;
time=toc;
t
end
if(Nb=='GL')%用复化Gauss_LegendreI
%先根据复化Gauss_LegendreI公式的余项估计步长
promps={'请输入用复化Gauss_LegendreI公式应取的步长:
'};
result=inputdlg(promps,'实验4.2',1,{'0.01'});
h=str2num(char(result));
if(h<=0)
errordlg('请输入正确的步长!
');
return;
end
tic;
N=floor((b-a)/h);t=0;
fork=0:
N-1
xk=a+k*h+h/2;
t=t+fun(xk-h/(2*sqrt(3)))+fun(xk+h/(2*sqrt(3)));
end
t=t*h/2;
time=toc;
t
end
switchNb_f
case1
disp('精确解:
ln2-ln3=-0.4054651081')
disp(['绝对误差:
',num2str(abs(t+0.4054651081))]);
disp(['运行时间:
',num2str(time)]);
case2
disp('精确解:
pi=3.979')
disp(['绝对误差:
',num2str(abs(t-pi))]);
disp(['运行时间:
',num2str(time)]);
case3
disp('精确解:
2/ln3=1.368')
disp(['绝对误差:
',num2str(abs(t-1.368))]);
disp(['运行时间:
',num2str(time)]);
case4
disp('精确解:
e^2=7.065')
disp(['绝对误差:
',num2str(abs(t-7.065))]);
disp(['运行时间:
',num2str(time)]);
end
1.当选用复化梯形公式时:
(1)式运行结果为:
t=-0.351
精确解:
ln2-ln3=-0.4054651081
绝对误差:
1.3944e-008
运行时间:
0.003
(2)式运行结果为:
t=3.336
精确解:
pi=3.979
绝对误差:
3.9736e-008
运行时间:
0.005
(3)式运行结果为:
t=1.861
精确解:
2/ln3=1.368
绝对误差:
4.3655e-008
运行时间:
0.016
(4)式运行结果为:
t=7.610
精确解:
e^2=7.065
绝对误差:
2.0775e-008
运行时间:
0.007
2.当选用复化Simpson公式进行计算时:
(1)式运行结果为:
t=-0.7519
精确解:
ln2-ln3=-0.4054651081
绝对误差:
2.7519e-011
运行时间:
0.022
(2)式运行结果为:
t=3.979
精确解:
pi=3.979
绝对误差:
0
运行时间:
0.021
(3)式运行结果为:
t=1.288
精确解:
2/ln3=1.368
绝对误差:
9.2018e-012
运行时间:
0.019
(4)式运行结果为:
t=7.118
精确解:
e^2=7.065
绝对误差:
9.0528e-011
运行时间:
0.021
3.当选用复化Gauss-LegendreI型公式进行计算时:
(1)式运行结果为:
t=-0.5262
精确解:
ln2-ln3=-0.4054651081
绝对误差:
4.7385e-012
运行时间:
0.023
(2)式运行结果为:
t=3.979
精确解:
pi=3.979
绝对误差:
1.3323e-015
运行时间:
0.021
(3)式运行结果为:
t=1.754
精确解:
2/ln3=1.368
绝对误差:
6.1431e-012
运行时间:
0.019
(4)式运行结果为:
t=7.046
精确解:
e^2=7.065
绝对误差:
1.441e-012
运行时间:
0.021
结果分析:
当选用复化梯形公式时,对步长的事前估计所要求的步长很小,选取的节点很多,误差绝对限要达到
时,对不同的函数n的取值需达到1000-10000之间,计算量是很大。
用复化simpson公式对步长的事前估计所要求的步长相对大些,选取的节点较少,误差绝对限要达到
时,对不同的函数n的取值只需在10-100之间,计算量相对小了很多,可满足用较少的节点达到较高的精度,比复化梯形公式的计算量小了很多。
用复化simpson公式计算所得的结果比用复化梯形公式计算所得的结果精度高很多,而且计算量小。
当选用Gauss-LagrangeI型公式进行计算时,选用较少的节点就可以达到很高的精度。
实验5.1常微分方程性态和R-K法稳定性试验
实验目的:
考察下面微分方程右端项中函数y前面的参数对方程性态的影响(它可使方程为好条件的或坏条件的)和研究计算步长对R-K法计算稳定性的影响。
实验容及要求:
实验题目:
常微分方程初值问题
其中,
。
其精确解为
。
实验要求:
本实验题都用4阶经典R-K法计算。
(1)对参数a分别取4个不同的数值:
一个大的正值,一个小的正值,一个绝对值小的负值和一个绝对值大的负值。
取步长h=0.01,分别用经典的R-K法计算,将四组计算结果画在同一图上,进行比较并说明相应初值问题的性态。
(2)取参数a为一个绝对值不大的负值和两个计算步长,一个步长使参数ah在经典R-K法的稳定域,另一个步长在经典R-K法的稳定域外。
分别用经典R-K法计算并比较计算结果。
取全域等距的10个点上的计算值,列表说明。
实验程序:
Matlab程序如下:
functioncharp5RK
%数值试验5.1:
常微分方程性态和R-K法稳定性试验
%输入:
参数a,步长h
%输出:
精确解和数值解图形对比
%clf;
result=inputdlg({'请输入[-50,50]间的参数a:
'},'实验5.1',1,{'-40'});
a=str2num(char(result));
if(a<-50|a>50)errordlg('请输入正确的参数a!
');return;end
result=inputdlg({'请输入(01)之间的步长:
'},'实验5.1',1,{'0.01'});
h=str2num(char(result));
if(h<=0|h>=1)errordlg('请输入正确的(01)间的步长!
');return;end
x=0:
h:
1;
y=x;
N=length(x);
y
(1)=1;
func=inline('1+(y-x).*a');
forn=1:
N-1
k1=func(a,x(n),y(n));
k2=func(a,x(n)+h/2,y(n)+k1*h/2);
k3=func(a,x(n)+h/2,y(n)+k2*h/2);
k4=func(a,x(n)+h,y(n)+k3*h);
y(n+1)=y(n)+h*(k1+2*k2+2*k3+k4)/6;