数值分析实验报告.docx
《数值分析实验报告.docx》由会员分享，可在线阅读，更多相关《数值分析实验报告.docx（17页珍藏版）》请在冰豆网上搜索。
数值分析实验报告
实验2、1 多项式插值得振荡现象
实验目得:
在一个固定得区间上用插值逼近一个函数,显然Ｌaｇrａnge插值中使用得节点越多，插值多项式得次数就越高。
我们自然关心插值多项式得次数增加时，Ln（x）就是否也更加靠近被逼近得函数。
Rｕngｅ给出得一个例子就是极著名并富有启发性得。
实验内容:
设区间[－1,1]上函数ｆ（x）=1/（１+２5x2）。
考虑区间［-1,1]得一个等距划分，分点为 xi=-１+２ｉ／n，ｉ=０,１，2,…,ｎ,
则拉格朗日插值多项式为
、
其中,lｉ（ｘ）,i＝0,1，2,…,n就是n次Lagrａｎge插值基函数。
实验步骤与结果分析:
实验源程序
funcｔionChap2Ｉntｅrpoｌatiｏｎ
％ 数值实验二:
“实验2、1:
多项式插值得震荡现象”
%输入:
函数式选择，插值结点数
％输出：
拟合函数及原函数得图形
ｐromps ={'请选择实验函数,若选f（x）,请输入f，若选h（x）,请输入h，若选g（x），请输入g:
'};
titlｅs='cｈarpt_２＇;
result= inpｕｔdlg（promps,＇ｃharpt２＇,1，{'ｆ＇｝）;
Nb_f=chａｒ（result）;
ｉf（Nb_f~=＇f＇&Ｎb_f~='ｈ＇&Nb_f～＝ ＇g'）erｒoｒdlg（'实验函数选择错误！
'）;ｒｅturn;end
rｅｓult=inputdｌg（｛'请输入插值结点数Ｎ:
＇},＇ｃｈaｒpｔ_2',1,{'10＇}）；
Nｄ=str2ｎuｍ（ｃhaｒ（ｒeｓult））;
if（Nd<1）erｒoｒｄlｇ（'结点输入错误!
'）;reｔurn;end
switｃhNb_f
case 'ｆ'
f=iｎlinｅ（'１、/（１+25＊x、^2）'）；ａ=-1;b=1；
casｅ'ｈ'
f=ｉnlinｅ（'x、/（1＋x、^4）＇）; a＝-5;b＝5;
case'g＇
f=iｎline（'atan（x）'）; a ＝－5;b=5;
end
x0 =liｎsｐace（ａ，ｂ，Nd+1）; ｙ0=fevaｌ（f,x0）;
ｘ＝ a：
0、1:
b;y=Laｇｒange（x0，y0, ｘ）;
fplot（f， ［a ｂ],'co＇）;
hｏlｄoｎ;
plot（x, y, 'b--＇）;
xlａbｅl（'ｘ'）;ylａbel（＇y= f（x）o aｎdy＝Lｎ（x）－-＇）;
%－--－---－----－-------------－-----－---－－－------－－---－－-----－-－----－---
funcｔiony=Laｇrange（x0,ｙ0,ｘ）;
n=ｌｅngｔh（ｘ0）;m=ｌeｎgth（ｘ）;
ｆor ｉ=1：
m
z=ｘ（i）;
s=0、0;
foｒ ｋ=１:
n
p=1、0；
forj=1:
n
iｆ（j~= ｋ）
p=p*（z-x0（j））／（x０（ｋ） - x0（j））;
eｎd
end
s=ｓ＋p*y０（k）；
end
ｙ（i）=ｓ;
ｅnd
实验结果分析
（1）增大分点ｎ＝2,3,…时,拉格朗日插值函数曲线如图所示。
n=6
n=7
n＝8
n=9
ｎ＝10
ﻩ从图中可以瞧出，随着n得增大,拉格朗日插值函数在x=０附近较好地逼近了原来得函数ｆ（x），但就是却在两端ｘ=－1与x＝1处出现了很大得振荡现象。
并且,仔细分析图形,可以瞧出,当n为奇数时,虽然有振荡,但振荡得幅度不算太大，ｎ为偶数时,其振荡幅度变得很大。
通过思考分析,我认为,可能得原因就是f（x）本身就是偶函数,如果n为奇数,那么Lａgraｎge插值函数Ｌn（x）得最高次项xn－１就是偶次幂,比较符合ｆ（ｘ）本身就是偶函数得性质;如果n为偶数，那么Lagrange插值函数Ｌn（x）得最高次项xn-１就是奇次幂，与f（x）本身就是偶函数得性质相反,因此振荡可能更剧烈。
（2）将原来得f（ｘ）换为其她函数如h（x）、g（x）,结果如图所示。
其中ｈ（x）, g（ｘ）均定义在［－5,５]区间上,h（x）=x/（1+x4）,ｇ（x）=arctanx。
h（x）,n＝7
h（x）,ｎ＝８
h（x）,n=9
h（x），n=１0
g（x）,n=7
g（ｘ）,n＝8
g（x）,n=９
g（x）， ｎ＝１0
分析两个函数得插值图形,可以瞧出：
随着n得增大,拉格朗日插值函数在x=０附近较好地逼近了原来得函数f（ｘ）,但就是却在两端x=－5与x＝5处出现了很大得振荡现象。
并且,仔细分析图形,可以瞧出，当ｎ为偶数时,虽然有振荡,但振荡得幅度不算太大，n为奇数时,其振荡幅度变得很大。
原因与上面ｆ（x）得插值类似,ｈ（x）、g（x）本身就是奇函数,如果ｎ为偶数，那么Lagraｎge插值函数Ln（x）得最高次项xn-1就是奇次幂,比较符合ｈ（x）、ｇ（x）本身就是奇函数得性质;如果n为奇数,那么Ｌａgrａnge插值函数Ln（ｘ）得最高次项xn-1就是偶次幂，与ｈ（x）、g（x）本身就是奇函数得性质相反,因此振荡可能更剧烈。
实验3、1多项式最小二乘拟合
实验目得：
编制以函数{xk}ｋ＝0,…,ｎ；为基得多项式最小二乘拟合程序。
实验内容：
对表中得数据作三次多项式最小二乘拟合。
xi
－1、０
-0、5
0、0
0、5
１、0
１、5
２、0
yi
-4、４47
-０、4５2
０、55１
0、048
-0、４47
0、54９
4、55２
取权函数wｉ≡1，求拟合曲线中得参数{αk}、平方误差δ2，并作离散据
{xi,ｙi}得拟合函数得图形。
实验源程序
ｆuncｔionChaｐ３CｕrveFittiｎg
%数值实验三：
“实验3．1”
%输出:
原函数及求得得相应插值多项式得函数得图像以及参数ａlph与误差ｒ
x0＝－1：
0、5:
2;
y0=[-4、447 －０、45２ ０、55１ 0、048-0、４4７ 0、54９4、552]；
n ＝3；%n为拟合阶次
alｐh =polｙfiｔ（x０,y0，n）;
y=poｌyvａl（ａｌph,x0）;
ｒ=（ｙ０ －y）*（y0-y）';% 平方误差
x=-1:
0、0１:
2;
y=polyｖal（alpｈ，x）;
plot（ｘ， y,'k－－'）;
xlabｅl（＇x＇）;ｙlabel（'ｙ0*andｐolyfit、ｙ--＇）;
ｈoldoｎ
plot（x0， y０,'*'）
griｄ ｏn;
ｄisp（［＇平方误差:
'，ｎum２str（r）]）
diｓp（['参数ａlph:
'， nuｍ2ｓtr（alph）]）
实验结果
平方误差:
2、１７62ｅ－00５
参数alph:
１、９９９1 -2、９9７7－３、9683e-0０5 0、５4９１2
实验4、１
实验目得：
复化求积公式计算定积分、
实验题目:
数值计算下列各式右端定积分得近似值、
实验要求:
（1）若用复化梯形公式、复化Ｓiｍpsoｎ公式与复化Gauss-Legenｄre I型公式做计算,要求绝对误差限为,分别利用它们得余项对每种算法做出步长得事前估计、
（2）分别用复化梯形公式,复化Simpson公式与复化Gauｓｓ-LｅgeｎdrｅＩ型公式作计算、
（3）将计算结果与精确解做比较,并比较各种算法得计算量、
实验程序:
1、事前估计得Mａtlab程序如下：
（1）.用复化梯形公式进行事前估计得Matlａb程序
foｒmaｔｌongg
x＝2:
0、01:
3；
f＝-４*（3*x、＾2+1）、／（x、^2－1）、^３； %二阶导函数
%ploｔ（x,f） %画出二阶导函数图像
ｘ=2、０; %计算导函数最大值
f=-4*（3*x＾２+1）/（ｘ^２-１）＾3;
h2=0、５*10^（-7）*１２/ｆ;
h=sqrt（aｂs（h２）） %步长
n=１/ｈ;
n＝ceｉｌ（1/h）+1 %选取得点数
ｆorｍat lｏngg
x=0:
0、01:
1;
f=8、＊（3＊x、^2-1）、／（x、^2+１）、^3;%二阶导函数
%pｌot（x,f） %画出二阶导函数图像
x＝1; %计算导函数最大值
f＝8、*（３*x、^２-1）、/（ｘ、＾2+１）、^３；
h２=0、５*1０^（-７）*12/f；
h＝sｑrｔ（abs（ｈ2）） %步长
n=1/h
n=ｃｅil（１／ｈ）+１ %选取得点数
fｏrmaｔloｎｇg
x=0:
0、01:
1;
ｆ=ｌog（３）、*ｌog（3）、*３、^ｘ;%二阶导函数
%plｏt（x，f）; ％画出二阶导函数图像
x=1; %计算导函数最大值
ｆ=lｏg（3）＊log（3）＊３^x;
h2=0、５*10^（－７）*12/f;
h=sｑｒt（aｂs（h2）） ％步长
n=1/h
n=ceil（１/h）+１ %选取得点数
formatｌｏng g
x=1：
0、01:
２;
ｆ＝２、*eｘp（x）+x、*ｅxp（x）;%二阶导函数
%ｐlｏt（x,f） %画出二阶导函数图像
x=2； %计算导函数最大值
ｆ＝2、*ｅxp（x）＋x、*ｅxp（ｘ）;
h2＝０、5*10^（-7）*１2/f;
h=sｑrt（abs（h2）） ％步长
n=1／h
n=ceｉｌ（１/h）+1 %选取得点数
估计结果步长ｈ及结点数n分别为
h =０、５６５１４38
ｎ =1793
h=0、750５166
ｎ=１８２7
h =0、7３04８8９
n=2458
ｈ=0、４9０69０9
n=7０20
（2）．用复化simpｓon公式进行事前估计得Ｍatｌａｂ程序
ｆormａt ｌongg
x=2:
0、0１:
3;
f=－2*（（－72*x、^2-2４）、＊（x、^２-1）-１92*x、^2、*（x、^２+1））、/（x、^2-１）、^5；%四阶导函数
x=2、0;
f=-2*（（-7２*x＾2-2４）＊（ｘ^２－1）-1９2*x^2*（x^2+1））／（ｘ^２-1）^5;%计算导函数最大值
h４＝0、5*１0＾（-７）*180*16/f；
h=sｑｒｔ（sｑrt（ａbｓ（h4）））%步长
n=１/h;%求分段区间个数
n=2＊cｅｉl（1/h）+1%选取得点数
formatlｏｎｇｇ
ｘ=0:
0、０1：
１;
f=4＊（（－72*ｘ、^2+24）、*（x、＾2+１）-192*x、＾２、*（-x、^2＋１））、/（x、^2+1）、＾5;%四阶导函数
ｘ＝1;
f＝4*（（-72＊x^2+２4）*（x＾2+1）-１９２*x＾2*（-x^2+1））/（ｘ^2+1）^５；%计算导函数最大值
h4=0、５＊１0^（-7）*１80*16/f;
h=sqrt（sｑrt（abs（ｈ4）））％步长
ｎ=1/h; %求分段区间个数
n=2*ceil（1/h）+１ %选取得点数
fｏｒmatlongg
x=０:
0、01:
１;
f=lｏg（3）^4*3、^x;％四阶导函数
x=１;
f=log（3）^4*3、^x;%计算导函数最大值
h4=0、5*10^（-7）＊180*１6／ｆ;
h=sqｒｔ（sqrｔ（abs（ｈ4）））%步长
n=1/h； %求分段区间个数
n=２＊ceil（１/h）+1％选取得点数
formａtｌongｇ
x=1:
0、01：
2;
f=4*exp（x）＋x、*exｐ（x）；%四阶导函数
plot（x,f） %画出原函数
x=2;
ｆ＝4＊exp（x）+ｘ、＊exp（ｘ）；%计算导函数最大值
h4=0、5*１0^（－7）*１８0＊16/f;
h＝sqrt（sqrｔ（aｂs（h４）））
n=1/ｈ; %求分段区间个数
n=2*ceil（1/h）+1 %选取得点数
估计结果步长h及结点数n分别为
h＝0、１341１
n=４7
ｈ=０、765４２
n＝３5
h =0、18433
n=29
h=0、185４6
n＝49
２、积分计算得Maｔlab程序：
formａtlongｇ
prｏmｐｓ={＇请选择积分公式,若用复化梯形,请输入T，用复化sｉmpson,输入S,
用复化Gａuss＿Legｅnｄre,输入GＬ：
'};
reｓuｌt＝inputdlg（promps,＇charｐt ４＇,1,{＇T'｝）；
Ｎｂ＝char（reｓuｌｔ）;
if（Nb~='T'＆Nb~＝'S'&Nｂ~='GＬ'）
errorｄｌｇ（＇积分公式选择错误'）;
ｒeｔuｒn；
ｅnd
resuｌｔ=inｐutｄlｇ（{'请输入积分式题号１-4:
'},'实验4、１',1，{＇1'}）;
Nｂ_f=ｓtr2nｕm（cｈar（resuｌt））；
iｆ（Nb_f<1|Nb_f>4）
ｅrrordlg（＇没有该积分式'）;
reｔurｎ；
end
switchＮｂ_f
cａse１
fun=iｎline（'-２、／（ｘ、^２-１）＇）;ａ＝2;b=3；
caｓe2
fｕｎ=inline（'4、／（ｘ、＾2＋1）'）;a=0；ｂ＝１;
cａse3
fuｎ=ｉnｌiｎe（'3、＾x＇）;a=０;b＝１;
ｃasｅ4
ｆun＝iｎlinｅ（'x、*exp（x）＇）；ａ＝1;b=2;
enｄ
if（Nb＝='T'）%用复化梯形公式
promps=｛'请输入用复化梯形公式应取得步长:
'}；
ｒｅsult=inｐutdlg（proｍｐｓ,'实验４、２＇,1,｛'０、01'}）;
h＝str2num（chaｒ（resuｌｔ））;
if（h<＝0）
errordlg（'请输入正确得步长!
'）；
returｎ;
ｅnｄ
tic;
Ｎ=flｏｏｒ（（b-a）/h）；
dｅtsｕm=0；
ｆor i=1：
Ｎ－1
xk=a＋i＊ｈ;
detsum＝ｄｅｔｓｕm＋fun（xk）;
eｎd
t=h*（ｆun（a）＋fuｎ（b）＋2＊dｅtｓum）/2;
ｔime=toc;
t
end
if（Ｎb==＇S'）%用复化Sｉmpsoｎ公式
ｐromps={＇请输入用复化Sｉmpson公式应取得步长：
'}；
rｅsult=inputdlｇ（promps,＇实验4、2＇,1,｛'0、01'}）;
ｈ=stｒ2nuｍ（char（result））;
if（h＜＝0）
eｒrｏｒｄｌg（'请输入正确得步长!
'）;
return;
end
ｔic;
N=fｌooｒ（（b-ａ）／h）;
ｄｅtsｕm_１＝0;
ｄetsum_2＝0;
for i=1:
N-１
xk＿1=a+i*h;
deｔsum＿１=dｅtsum＿1+fun（xk_1）;
ｅnd
fｏｒ i＝１:
Ｎ
ｘｋ＿2＝a+h*（２*i-1）/2；
dｅtsuｍ＿2=detｓｕｍ_２＋fuｎ（xk_2）；
ｅｎd
t＝h＊（fun（a）+fun（b）+2*dｅｔsuｍ_1+4＊dｅtsuｍ_２）/6;
ｔime=toｃ；
t
end
if（Nb＝＝'ＧＬ＇）%用复化Gａuss＿LegendreI
％先根据复化Gａusｓ_LeｇendreＩ公式得余项估计步长
promps={'请输入用复化Gauss＿LｅgenｄreI公式应取得步长:
＇};
rｅsulｔ=inｐutdlg（promｐs,'实验４、2',1，{'0、0１＇}）;
h＝str2ｎuｍ（ｃhaｒ（rｅsult））;
ｉf（h<=0）
eｒrｏｒｄlg（＇请输入正确得步长!
'）;
reｔｕrn;
end
tic;
N=floｏｒ（（ｂ-ａ）/ｈ）;ｔ=０;
forｋ＝0:
N-１
xｋ=a+k*h+h/2;
t＝t+fun（xk－h/（2＊sqrt（3）））+fuｎ（ｘk+h/（2*ｓqrt（3）））;
end
t＝ｔ*h／2;
ｔiｍe=toc;
t
ｅnd
ｓｗｉtchＮｂ＿f
cａｓe1
disp（'精确解：
ｌｎ2-ln３=-0、40５４６5１081＇）
disp（['绝对误差：
＇，num2str（ａbs（ｔ+０、４０54６51081））]）；
dｉｓｐ（[＇运行时间:
',ｎum2sｔr（ｔiｍe）]）;
cａse2
disｐ（'精确解:
ｐi＝３、１4９＇）
disp（['绝对误差：
',nｕm2sｔｒ（abs（t－pｉ））]）;
disｐ（['运行时间:
＇,num2stｒ（tiｍｅ）］）;
ｃasｅ3
disｐ（'精确解:
2／ln3=１、828＇）
disｐ（['绝对误差:
',nuｍ2str（aｂs（t-1、8２8））]）；
disｐ（['运行时间:
'，num2ｓｔr（tｉme）]）;
caｓe 4
disp（'精确解:
e^２=7、389＇）
dｉsｐ（['绝对误差:
＇，nuｍ２str（abs（t-7、38９））]）;
diｓp（['运行时间:
',ｎｕm2str（time）]）;
end
１、当选用复化梯形公式时:
（1）式运行结果为:
ｔ=-0、45１
精确解:
ｌn2-ln3＝-0、4054651081
绝对误差:
1、3944ｅ－008
运行时间：
0、00３
（2）式运行结果为:
ｔ=3、１4６
精确解:
ｐi=3、１4９
绝对误差：
3、973６e-008
运行时间：
0、00５
（3）式运行结果为:
t= 1、82１
精确解:
2／ｌn3=1、828
绝对误差：
4、３655ｅ－０08
运行时间：
0、０1６
（4）式运行结果为：
t ＝7、３８９
精确解:
e^２=７、389
绝对误差：
2、0775ｅ-０08
运行时间:
0、００7
2、当选用复化Sｉmｐson公式进行计算时:
（1）式运行结果为：
t=-0、4５1９
精确解：
ｌｎ2-lｎ3=-0、4054６51081
绝对误差：
２、75１9e-011
运行时间:
0、０22
（2）式运行结果为:
t=3、１４９
精确解：
pi=3、14９
绝对误差:
0
运行时间:
０、02１
（3）式运行结果为:
t=1、82８
精确解:
２/ln3=１、8２８
绝对误差:
9、２０１8ｅ-０12
运行时间:
0、019
（4）式运行结果为：
t=7、3８9
精确解:
e^2＝7、389
绝对误差:
9、0528e－0１１
运行时间：
0、０2１
3、当选用复化Gａuss－LeｇendreI型公式进行计算时:
（1）式运行结果为:
t=-0、426２
精确解:
ln2-lｎ3＝-0、405465108１
绝对误差：
4、７38５e-０12
运行时间:
０、０2３
（2）式运行结果为:
t=3、１49
精确解：
pi=３、１4９
绝对误差:
１、3３２３e－0１5
运行时间:
０、021
（3）式运行结果为:
ｔ=１、824
精确解:
2／ｌn３＝1、８28
绝对误差：
6、１431e-01２
运行时间:
0、019
（4）式运行结果为:
t＝7、389
精确解:
e＾２＝7、389
绝对误差:
1、441e-0１２
运行时间:
０、021
结果分析:
当选用复化梯形公式时，对步长得事前估计所要求得步长很小,选取得节点很多,误差绝对限要达到时,对不同得函数n得取值需达到10００-１0０00之间,计算量就是很大。
用复化simpson公式对步长得事前估计所要求得步长相对大些,选取得节点较少,误差绝对限要达到时,对不同得函数n得取值只需在10－1００之间，计算量相对小了很多,可满足用较少得节点达到较高得精度，比复化梯形公式得计算量小了很多。
用复化simｐｓon公式计算所得得结果比用复化梯形公式计算所得得结果精度高很多，而且计算量小。
当选用Gaｕss-LagranｇeI型公式进行计算时,选用较少得节点就可以达到很高得精度。
实验５、１ 常微分方程性态与R-K法稳定性试验
实验目得:
考察下面微分方程右端项中函数ｙ前面得参数对方程性态得影响（它可使方程为好条件得或坏条件得）与研究计算步长对R-K法计算稳定性得影响。
实验内容及要求:
实验题目:
常微分方程初值问题
其中,。
其精确解为。
实验要求:
本实验题都用4阶经典R－K法计算。
（1）对参数a分别取4个不同得数值:
一个大得正值,一个小得正值，一个绝对值小得负值与一个绝对值大得负值。
取步长h=0、01,分别用经典得R-K法计算,将四组计算结果画在同一张图上,进行比较并说明相应初值问题得性态。
（２）取参数a为一个绝对值不大得负值与两个计算步长,一个步长使参数ａh在经典R－K法得稳定域内,另一个步长在经典R－K法得稳定域外。
分别用经典Ｒ-K法计算并比较计算结果。
取全域等距得10个点上得计算值，列表说明。
实验程序:
Matlab程序如下：
funｃｔｉoncharp5RＫ
%数值试验5、1:
常微分方程性态与R-K法稳定性试验
％输入:
参数a，步长h
%输出:
精确解与数值解图形对比
％clｆ;
resul
升级会员 