数值分析实验报告.docx
《数值分析实验报告.docx》由会员分享,可在线阅读,更多相关《数值分析实验报告.docx(17页珍藏版)》请在冰豆网上搜索。
数值分析实验报告
实验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, [a b],'co');
holdon;
plot(x, y, 'b--');
xlabel('x');ylabel('y= f(x)o andy=Ln(x)--');
%--------------------------------------------------------------------
functiony=Lagrange(x0,y0,x);
n=length(x0);m=length(x);
for i=1:
m
z=x(i);
s=0、0;
for k=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、452 0、551 0、048-0、447 0、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,'*')
grid on;
disp(['平方误差:
',num2str(r)])
disp(['参数alph:
', num2str(alph)])
实验结果
平方误差:
2、1762e-005
参数alph:
1、9991 -2、9977-3、9683e-005 0、54912
实验4、1
实验目得:
复化求积公式计算定积分、
实验题目:
数值计算下列各式右端定积分得近似值、
实验要求:
(1)若用复化梯形公式、复化Simpson公式与复化Gauss-Legendre I型公式做计算,要求绝对误差限为,分别利用它们得余项对每种算法做出步长得事前估计、
(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 %选取得点数
format longg
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 %选取得点数
formatlong g
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、5651438
n =1793
h=0、7505166
n=1827
h =0、7304889
n=2458
h=0、4906909
n=7020
(2).用复化simpson公式进行事前估计得Matlab程序
format longg
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,'charpt 4',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;
for i=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;
for i=1:
N-1
xk_1=a+i*h;
detsum_1=detsum_1+fun(xk_1);
end
for i=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、149')
disp(['绝对误差:
',num2str(abs(t-pi))]);
disp(['运行时间:
',num2str(time)]);
case3
disp('精确解:
2/ln3=1、828')
disp(['绝对误差:
',num2str(abs(t-1、828))]);
disp(['运行时间:
',num2str(time)]);
case 4
disp('精确解:
e^2=7、389')
disp(['绝对误差:
',num2str(abs(t-7、389))]);
disp(['运行时间:
',num2str(time)]);
end
1、当选用复化梯形公式时:
(1)式运行结果为:
t=-0、451
精确解:
ln2-ln3=-0、4054651081
绝对误差:
1、3944e-008
运行时间:
0、003
(2)式运行结果为:
t=3、146
精确解:
pi=3、149
绝对误差:
3、9736e-008
运行时间:
0、005
(3)式运行结果为:
t= 1、821
精确解:
2/ln3=1、828
绝对误差:
4、3655e-008
运行时间:
0、016
(4)式运行结果为:
t =7、389
精确解:
e^2=7、389
绝对误差:
2、0775e-008
运行时间:
0、007
2、当选用复化Simpson公式进行计算时:
(1)式运行结果为:
t=-0、4519
精确解:
ln2-ln3=-0、4054651081
绝对误差:
2、7519e-011
运行时间:
0、022
(2)式运行结果为:
t=3、149
精确解:
pi=3、149
绝对误差:
0
运行时间:
0、021
(3)式运行结果为:
t=1、828
精确解:
2/ln3=1、828
绝对误差:
9、2018e-012
运行时间:
0、019
(4)式运行结果为:
t=7、389
精确解:
e^2=7、389
绝对误差:
9、0528e-011
运行时间:
0、021
3、当选用复化Gauss-LegendreI型公式进行计算时:
(1)式运行结果为:
t=-0、4262
精确解:
ln2-ln3=-0、4054651081
绝对误差:
4、7385e-012
运行时间:
0、023
(2)式运行结果为:
t=3、149
精确解:
pi=3、149
绝对误差:
1、3323e-015
运行时间:
0、021
(3)式运行结果为:
t=1、824
精确解:
2/ln3=1、828
绝对误差:
6、1431e-012
运行时间:
0、019
(4)式运行结果为:
t=7、389
精确解:
e^2=7、389
绝对误差:
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;
resul