數(shù)值分析實(shí)驗(yàn)報(bào)告

1、實(shí)驗(yàn) 2.1 多項(xiàng)式插值的振蕩現(xiàn)象實(shí)驗(yàn)?zāi)康模涸谝粋€(gè)固定的區(qū)間上用插值逼近一個(gè)函數(shù)，顯然Lagrange 插值中使用的節(jié)點(diǎn)越多，插值多項(xiàng)式的次數(shù)就越高。我們自然關(guān)心插值多項(xiàng)式的次數(shù)增加時(shí)，Ln(x) 是否也更加靠近被逼近的函數(shù)。 Runge 給出的一個(gè)例子是極著名并富有啟發(fā)性的。實(shí)驗(yàn)內(nèi)容：設(shè)區(qū)間 -1， 1上函數(shù) f(x)=1/(1+25x2)?？紤]區(qū)間 -1,1 的一個(gè)等距劃分，分點(diǎn)為xi= -1 + 2i/n ， i=0 ， 1， 2， , ， n，則拉格朗日插值多項(xiàng)式為n1Ln (x)li ( x) .i 0 1 25xi2其中， li(x) ， i=0 ，1， 2， , ， n 是 n

2、次 Lagrange 插值基函數(shù)。實(shí)驗(yàn)步驟與結(jié)果分析：實(shí)驗(yàn)源程序function Chap2Interpolation% 數(shù)值實(shí)驗(yàn)二： “實(shí)驗(yàn) 2.1：多項(xiàng)式插值的震蕩現(xiàn)象”% 輸入：函數(shù)式選擇，插值結(jié)點(diǎn)數(shù)% 輸出：擬合函數(shù)及原函數(shù)的圖形promps = ' 請(qǐng)選擇實(shí)驗(yàn)函數(shù)，若選f(x), 請(qǐng)輸入 f, 若選 h(x), 請(qǐng)輸入 h,若選 g(x), 請(qǐng)輸入titles = 'charpt_2'result = inputdlg(promps,'charpt 2',1,'f');Nb_f = char(result);if(Nb_f =

3、'f' & Nb_f = 'h' & Nb_f = 'g')errordlg('實(shí)驗(yàn)函數(shù)選擇錯(cuò)誤！');return;endresult = inputdlg(' 請(qǐng)輸入插值結(jié)點(diǎn)數(shù)N:','charpt_2',1,'10');g:'Nd = str2num(char(result);if(Nd <1)errordlg(' 結(jié)點(diǎn)輸入錯(cuò)誤！');return;endswitch Nb_fcase 'f'f=inline(

4、9;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;endx0 = linspace(a, b, Nd+1); y0 = feval(f, x0);x = a:0.1:b; y = Lagrange(x0, y0, x);fplot(f, a b, 'co');hold on;plot(x, y, 'b-'

5、);xlabel('x'); ylabel('y = f(x) o and y = Ln(x)-');%-function y=Lagrange(x0, y0, x);n= length(x0); m=length(x);for i=1:mz=x(i);s=0.0;for k=1:np=1.0;for j=1:nif(j = k)p = p*(z - x0(j)/(x0(k) - x0(j);endends = s + p*y0(k);endy(i) = s;end實(shí)驗(yàn)結(jié)果分析(1) 增大分點(diǎn) n=2，3， , 時(shí)，拉格朗日插值函數(shù)曲線(xiàn)如圖所示。 n=6n=7n

6、=8n=9n=10從圖中可以看出，隨著n 的增大，拉格朗日插值函數(shù)在x=0 附近較好地逼近了原來(lái)的函數(shù) f(x) ，但是卻在兩端x= -1 和 x=1 處出現(xiàn)了很大的振蕩現(xiàn)象。并且，仔細(xì)分析圖形， 可以看出， 當(dāng) n 為奇數(shù)時(shí)， 雖然有振蕩， 但振蕩的幅度不算太大，n 為偶數(shù)時(shí)，其振蕩幅度變得很大。通過(guò)思考分析，我認(rèn)為，可能的原因是f(x) 本身是偶函數(shù)，如果 n 為奇數(shù)，那么Lagrange 插值函數(shù) L n(x) 的最高次項(xiàng)xn-1 是偶次冪，比較符合 f(x)本身是偶函數(shù)的性質(zhì)；如果n 為偶數(shù)，那么 Lagrange 插值函數(shù)L n(x) 的最高次項(xiàng) xn-1 是奇次冪，與 f(x) 本

7、身是偶函數(shù)的性質(zhì)相反，因此振蕩可能更劇烈。(2) 將原來(lái)的f(x) 換為其他函數(shù)如h(x)、 g(x) ，結(jié)果如圖所示。其中 h(x), g(x) 均定義在 -5， 5區(qū)間上， h(x)=x/(1+x 4)， g(x)=arctan x 。h(x), n=7h(x), n=8h(x), n=9h(x), n=10g(x), n=7g(x), n=8g(x), n=9g(x), n=10分析兩個(gè)函數(shù)的插值圖形，可以看出：隨著 n 的增大， 拉格朗日插值函數(shù)在x=0附近較好地逼近了原來(lái)的函數(shù)f(x) ，但是卻在兩端x= -5和x=5處出現(xiàn)了很大的振蕩現(xiàn)象。并且，仔細(xì)分析圖形， 可以看出， 當(dāng) n

8、為偶數(shù)時(shí)， 雖然有振蕩， 但振蕩的幅度不算太大，n 為奇數(shù)時(shí)，其振蕩幅度變得很大。原因和上面f(x) 的插值類(lèi)似，h(x) 、 g(x) 本身是奇函數(shù)，如果 n 為偶數(shù)， 那么 Lagrange 插值函數(shù) L n(x) 的最高次項(xiàng) xn-1 是奇次冪， 比較符合 h(x) 、g(x)本身是奇函數(shù)的性質(zhì)；如果 n 為奇數(shù)，那么 Lagrange 插值函數(shù) L n(x) 的最高次項(xiàng) xn-1 是偶次冪，與 h(x) 、 g(x) 本身是奇函數(shù)的性質(zhì)相反，因此振蕩可能更劇烈。實(shí)驗(yàn) 3.1 多項(xiàng)式最小二乘擬合實(shí)驗(yàn)?zāi)康模壕幹埔院瘮?shù) xkk=0, ,n；為基的多項(xiàng)式最小二乘擬合程序。實(shí)驗(yàn)內(nèi)容：對(duì)表中的數(shù)據(jù)

9、作三次多項(xiàng)式最小二乘擬合。xi-1.0-0.50.00.51.01.52.0yi-4.447-0.4520.5510.048-0.4470.5494.552n取權(quán)函數(shù) w i 1，求擬合曲線(xiàn)*k* xk中的參數(shù) k 、平方誤差 2，并作離散據(jù)k 0x i , yi 的擬合函數(shù)的圖形。實(shí)驗(yàn)源程序function Chap3CurveFitting% 數(shù)值實(shí)驗(yàn)三 :“實(shí)驗(yàn) 3.1”%輸出 :原函數(shù)及求得的相應(yīng)插值多項(xiàng)式的函數(shù)的圖像以及參數(shù)alph 和誤差 rx0= -1:0.5:2;y0= -4.447 -0.452 0.551 0.048 -0.447 0.549 4.552;n = 3; %

10、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 * and polyfit.y-');hold onplot(x0, y0, '*')grid on;disp(' 平方誤差 :', num2str(r)disp(' 參數(shù) alph:

11、9;, num2str(alph)實(shí)驗(yàn)結(jié)果平方誤差 :2.1762e-005參數(shù) alph:1.9991-2.9977 -3.9683e-0050.54912實(shí)驗(yàn) 4.1實(shí)驗(yàn)?zāi)康模簭?fù)化求積公式計(jì)算定積分.實(shí)驗(yàn)題目：數(shù)值計(jì)算下列各式右端定積分的近似值.實(shí)驗(yàn)要求：（1）若用復(fù)化梯形公式、復(fù)化Simpson 公式和復(fù)化Gauss-Legendre I 型公式做計(jì)算，要求1*107絕對(duì)誤差限為2，分別利用它們的余項(xiàng)對(duì)每種算法做出步長(zhǎng)的事前估計(jì).（ 2）分別用復(fù)化梯形公式，復(fù)化Simpson 公式和復(fù)化 Gauss-Legendre I 型公式作計(jì)算 .（3）將計(jì)算結(jié)果與精確解做比較，并比較各種算法的計(jì)

12、算量.實(shí)驗(yàn)程序：1.事前估計(jì)的Matlab 程序如下：（1）用復(fù)化梯形公式進(jìn)行事前估計(jì)的Matlab 程序format long gx=2:0.01:3;f=-4*(3*x.2+1)./(x.2-1).3; %二階導(dǎo)函數(shù)%plot(x,f)%畫(huà)出二階導(dǎo)函數(shù)圖像x=2.0;%計(jì)算導(dǎo)函數(shù)最大值f=-4*(3*x2+1)/(x2-1)3;h2=0.5*10(-7)*12/f;h=sqrt(abs(h2)%步長(zhǎng)n=1/h;n=ceil(1/h)+1%選取的點(diǎn)數(shù)format long gx=0:0.01:1;f=8.*(3*x.2-1)./(x.2+1).3;%二階導(dǎo)函數(shù)%plot(x,f)%畫(huà)出二階導(dǎo)

13、函數(shù)圖像x=1;% 計(jì)算導(dǎo)函數(shù)最大值f=8.*(3*x.2-1)./(x.2+1).3;h2=0.5*10(-7)*12/f;h=sqrt(abs(h2)%步長(zhǎng)n=1/hn=ceil(1/h)+1%選取的點(diǎn)數(shù)format long gx=0:0.01:1;f=log(3).*log(3).*3.x;%plot(x,f);x=1;% 二階導(dǎo)函數(shù)% 畫(huà)出二階導(dǎo)函數(shù)圖像% 計(jì)算導(dǎo)函數(shù)最大值f=log(3)*log(3)*3x;h2=0.5*10(-7)*12/f;h=sqrt(abs(h2)%步長(zhǎng)n=1/hn=ceil(1/h)+1% 選取的點(diǎn)數(shù)format long gx=1:0.01:2;f=2

14、.*exp(x)+x.*exp(x);%二階導(dǎo)函數(shù)%plot(x,f)x=2;%畫(huà)出二階導(dǎo)函數(shù)圖像% 計(jì)算導(dǎo)函數(shù)最大值f=2.*exp(x)+x.*exp(x);h2=0.5*10(-7)*12/f;h=sqrt(abs(h2)%步長(zhǎng)n=1/hn=ceil(1/h)+1%選取的點(diǎn)數(shù)估計(jì)結(jié)果步長(zhǎng)h 及結(jié)點(diǎn)數(shù)n 分別為h = 0.000558156305651438n =1793h =0.000547722557505166 n =1827h =0.000407071357304889 n =2458h =0.000142479094906909n =7020（2）用復(fù)化 simpson 公式進(jìn)行

15、事前估計(jì)的Matlab 程序format long gx=2:0.01:3;f=-2*(-72*x.2-24).*(x.2-1)-192*x.2.*(x.2+1)./(x.2-1).5;%四階導(dǎo)函數(shù)x=2.0;f=-2*(-72*x2-24)*(x2-1)-192*x2*(x2+1)/(x2-1)5; %計(jì)算導(dǎo)函數(shù)最大值h4=0.5*10(-7)*180*16/f;h=sqrt(sqrt(abs(h4)%步長(zhǎng)n=1/h;%求分段區(qū)間個(gè)數(shù)n=2*ceil(1/h)+1 % 選取的點(diǎn)數(shù)format long gx=0:0.01:1;f=4*(-72*x.2+24).*(x.2+1)-192*x.2

16、.*(-x.2+1)./(x.2+1).5;%四階導(dǎo)函數(shù)x=1;f=4*(-72*x2+24)*(x2+1)-192*x2*(-x2+1)/(x2+1)5; %計(jì)算導(dǎo)函數(shù)最大值h4=0.5*10(-7)*180*16/f;h=sqrt(sqrt(abs(h4)% 步長(zhǎng)n=1/h;% 求分段區(qū)間個(gè)數(shù)n=2*ceil(1/h)+1 % 選取的點(diǎn)數(shù)format long gx=0:0.01:1;f=log(3)4*3.x;%四階導(dǎo)函數(shù)x=1;f=log(3)4*3.x;%計(jì)算導(dǎo)函數(shù)最大值h4=0.5*10(-7)*180*16/f;h=sqrt(sqrt(abs(h4)% 步長(zhǎng)n=1/h;% 求分段

17、區(qū)間個(gè)數(shù)n=2*ceil(1/h)+1 % 選取的點(diǎn)數(shù)format long gx=1:0.01:2;f=4*exp(x)+x.*exp(x);%四階導(dǎo)函數(shù)plot(x,f) % 畫(huà)出原函數(shù)x=2;f=4*exp(x)+x.*exp(x); %計(jì)算導(dǎo)函數(shù)最大值h4=0.5*10(-7)*180*16/f;h=sqrt(sqrt(abs(h4)n=1/h;% 求分段區(qū)間個(gè)數(shù)n=2*ceil(1/h)+1 % 選取的點(diǎn)數(shù)估計(jì)結(jié)果步長(zhǎng)h 及結(jié)點(diǎn)數(shù)n 分別為h =0.0437490486013411 n =47h =0.0588566191276542 n =35h =0.07576451662184

18、33 n =29h =0.0424527247118546 n =492. 積分計(jì)算的 Matlab 程序：format long gpromps=' 請(qǐng)選擇積分公式，若用復(fù)化梯形，請(qǐng)輸入 T ，用復(fù)化 simpson，輸入 S，用復(fù)化 Gauss_Legendre，輸入 GL ： 'result=inputdlg(promps,'charpt 4',1,'T');Nb=char(result);if(Nb='T'&Nb='S'&Nb='GL')errordlg(' 積分公式

19、選擇錯(cuò)誤');return;endresult=inputdlg(' 請(qǐng)輸入積分式題號(hào)1-4： ',' 實(shí)驗(yàn) 4.1',1,'1');Nb_f=str2num(char(result);if(Nb_f<1|Nb_f>4)errordlg(' 沒(méi)有該積分式 ');return;endswitch Nb_fcase 1fun=inline('-2./(x.2-1)');a=2;b=3;case 2fun=inline('4./(x.2+1)');a=0;b=1;case 3fun=i

20、nline('3.x');a=0;b=1;case 4fun=inline('x.*exp(x)');a=1;b=2;endif(Nb='T')% 用復(fù)化梯形公式promps=' 請(qǐng)輸入用復(fù)化梯形公式應(yīng)取的步長(zhǎng)：'result=inputdlg(promps,' 實(shí)驗(yàn) 4.2',1,'0.01');h=str2num(char(result);if(h<=0)errordlg(' 請(qǐng)輸入正確的步長(zhǎng)！');return;endtic;N=floor(b-a)/h);detsum=

21、0;for i=1:N-1xk=a+i*h;detsum=detsum+fun(xk);endt=h*(fun(a)+fun(b)+2*detsum)/2;time=toc;tendif(Nb='S')% 用復(fù)化 Simpson 公式promps=' 請(qǐng)輸入用復(fù)化Simpson 公式應(yīng)取的步長(zhǎng)：'result=inputdlg(promps,'實(shí)驗(yàn)4.2',1,'0.01');h=str2num(char(result);if(h<=0)errordlg(' 請(qǐng)輸入正確的步長(zhǎng)！');return;endtic

22、;N=floor(b-a)/h);detsum_1=0;detsum_2=0;for i=1:N-1xk_1=a+i*h;detsum_1=detsum_1+fun(xk_1);endfor i=1:Nxk_2=a+h*(2*i-1)/2;detsum_2=detsum_2+fun(xk_2);endt=h*(fun(a)+fun(b)+2*detsum_1+4*detsum_2)/6;time=toc;tendif(Nb='GL')% 用復(fù)化 Gauss_Legendre I%先根據(jù)復(fù)化Gauss_Legendre I 公式的余項(xiàng)估計(jì)步長(zhǎng)promps=' 請(qǐng)輸入用復(fù)化

23、 Gauss_Legendre I 公式應(yīng)取的步長(zhǎng)： ' result=inputdlg(promps,' 實(shí)驗(yàn) 4.2',1,'0.01'); h=str2num(char(result);if(h<=0)errordlg(' 請(qǐng)輸入正確的步長(zhǎng)!');return;endtic;N=floor(b-a)/h);t=0;for k=0:N-1xk=a+k*h+h/2;t=t+fun(xk-h/(2*sqrt(3)+fun(xk+h/(2*sqrt(3);endt=t*h/2;time=toc;tendswitch Nb_fcase

24、1disp('精確解： ln2-ln3=-0.4054651081')disp(' 絕對(duì)誤差： ',num2str(abs(t+0.4054651081);disp(' 運(yùn)行時(shí)間： ',num2str(time);case 2disp('精確解： pi=3.14159265358979')disp(' 絕對(duì)誤差： ',num2str(abs(t-pi);disp(' 運(yùn)行時(shí)間： ',num2str(time);case 3disp('精確解： 2/ln3=1.82047845325368&#

25、39;)disp(' 絕對(duì)誤差： ',num2str(abs(t-1.82047845325368);disp(' 運(yùn)行時(shí)間： ',num2str(time);case 4disp('精確解： e2=7.38905609893065')disp(' 絕對(duì)誤差： ',num2str(abs(t-7.38905609893065);disp(' 運(yùn)行時(shí)間： ',num2str(time);end1. 當(dāng)選用復(fù)化梯形公式時(shí)：（1）式運(yùn)行結(jié)果為：t =-0.40546512204351精確解： ln2-ln3=-0.405

26、4651081絕對(duì)誤差： 1.3944e-008運(yùn)行時(shí)間： 0.003（2）式運(yùn)行結(jié)果為：t =3.14159261385336精確解： pi=3.14159265358979絕對(duì)誤差： 3.9736e-008運(yùn)行時(shí)間： 0.005（3）式運(yùn)行結(jié)果為：t = 1.82047849690861精確解： 2/ln3=1.82047845325368絕對(duì)誤差： 4.3655e-008運(yùn)行時(shí)間： 0.016（4）式運(yùn)行結(jié)果為：t =7.38905611970610精確解： e2=7.38905609893065絕對(duì)誤差： 2.0775e-008運(yùn)行時(shí)間： 0.0072. 當(dāng)選用復(fù)化 Simpson 公

27、式進(jìn)行計(jì)算時(shí)：（1）式運(yùn)行結(jié)果為：t =-0.405465108127519精確解： ln2-ln3=-0.4054651081絕對(duì)誤差： 2.7519e-011運(yùn)行時(shí)間： 0.022（2）式運(yùn)行結(jié)果為：t =3.14159265358979精確解： pi=3.14159265358979絕對(duì)誤差： 0運(yùn)行時(shí)間： 0.021（3）式運(yùn)行結(jié)果為：t =1.82047845326288精確解： 2/ln3=1.82047845325368絕對(duì)誤差： 9.2018e-012運(yùn)行時(shí)間： 0.019（4）式運(yùn)行結(jié)果為：t =7.38905609902118精確解： e2=7.38905609893065

28、絕對(duì)誤差： 9.0528e-011運(yùn)行時(shí)間： 0.0213. 當(dāng)選用復(fù)化 Gauss-Legendre I 型公式進(jìn)行計(jì)算時(shí)：（1）式運(yùn)行結(jié)果為：t =-0.405465108095262精確解： ln2-ln3=-0.4054651081絕對(duì)誤差： 4.7385e-012運(yùn)行時(shí)間： 0.023（2）式運(yùn)行結(jié)果為：t =3.14159265358979精確解： pi=3.14159265358979絕對(duì)誤差： 1.3323e-015運(yùn)行時(shí)間： 0.021（3）式運(yùn)行結(jié)果為：t =1.82047845324754精確解： 2/ln3=1.82047845325368絕對(duì)誤差： 6.1431e-0

29、12運(yùn)行時(shí)間： 0.019（4）式運(yùn)行結(jié)果為：t =7.38905607315046精確解： e2=7.38905609893065絕對(duì)誤差： 1.441e-012運(yùn)行時(shí)間： 0.021結(jié)果分析：當(dāng)選用復(fù)化梯形公式時(shí)， 對(duì)步長(zhǎng)的事前估計(jì)所要求的步長(zhǎng)很小，選取的節(jié)點(diǎn)很多， 誤差絕對(duì)限要達(dá)到1*1072時(shí)，對(duì)不同的函數(shù) n 的取值需達(dá)到 1000-10000 之間，計(jì)算量是很大。用復(fù)化 simpson 公式對(duì)步長(zhǎng)的事前估計(jì)所要求的步長(zhǎng)相對(duì)大些，選取的節(jié)點(diǎn)較少， 誤差絕對(duì)17限要達(dá)到* 10時(shí)，對(duì)不同的函數(shù) n 的取值只需在 10-100之間，計(jì)算量相對(duì)小了很多，2可滿(mǎn)足用較少的節(jié)點(diǎn)達(dá)到較高的精度，

30、比復(fù)化梯形公式的計(jì)算量小了很多。用復(fù)化simpson公式計(jì)算所得的結(jié)果比用復(fù)化梯形公式計(jì)算所得的結(jié)果精度高很多，而且計(jì)算量小。當(dāng)選用 Gauss-Lagrange I 型公式進(jìn)行計(jì)算時(shí)，選用較少的節(jié)點(diǎn)就可以達(dá)到很高的精度。實(shí)驗(yàn) 5.1常微分方程性態(tài)和 R-K 法穩(wěn)定性試驗(yàn)實(shí)驗(yàn)?zāi)康模嚎疾煜旅嫖⒎址匠逃叶隧?xiàng)中函數(shù)y 前面的參數(shù)對(duì)方程性態(tài)的影響 （它可使方程為好條件的或壞條件的）和研究計(jì)算步長(zhǎng)對(duì)R-K 法計(jì)算穩(wěn)定性的影響。實(shí)驗(yàn)內(nèi)容及要求：實(shí)驗(yàn)題目：常微分方程初值問(wèn)題'yayax1,0x1,其中，50a50 。其精確解為 y( x)eaxx 。實(shí)驗(yàn)要求：本實(shí)驗(yàn)題都用4 階經(jīng)典 R-K 法計(jì)算。

31、（ 1）對(duì)參數(shù) a 分別取 4 個(gè)不同的數(shù)值：一個(gè)大的正值，一個(gè)小的正值，一個(gè)絕對(duì)值小的負(fù)值和一個(gè)絕對(duì)值大的負(fù)值。取步長(zhǎng) h=0.01，分別用經(jīng)典的 R-K 法計(jì)算，將四組計(jì)算結(jié)果畫(huà)在同一張圖上，進(jìn)行比較并說(shuō)明相應(yīng)初值問(wèn)題的性態(tài)。（ 2）取參數(shù) a 為一個(gè)絕對(duì)值不大的負(fù)值和兩個(gè)計(jì)算步長(zhǎng)， 一個(gè)步長(zhǎng)使參數(shù) ah 在經(jīng)典 R-K 法的穩(wěn)定域內(nèi)，另一個(gè)步長(zhǎng)在經(jīng)典 R-K 法的穩(wěn)定域外。分別用經(jīng)典 R-K法計(jì)算并比較計(jì)算結(jié)果。取全域等距的10 個(gè)點(diǎn)上的計(jì)算值，列表說(shuō)明。實(shí)驗(yàn)程序：Matlab 程序如下：function charp5RK%數(shù)值試驗(yàn)5.1：常微分方程性態(tài)和R-K 法穩(wěn)定性試驗(yàn)%輸入：參數(shù)a，步長(zhǎng) h%輸出：精確解和數(shù)值解圖形對(duì)比%clf;result=inputdlg(' 請(qǐng)輸入 -50 ， 50間的參數(shù)a:',' 實(shí)驗(yàn) 5.1',1,'-40');a=str2num(char(result);if (a<-50|a>50)errordlg(' 請(qǐng)輸入正確的參數(shù)a!'); return;endresult=inputdlg(' 請(qǐng)輸入（ 0 1）之間的步長(zhǎng):',' 實(shí)驗(yàn) 5.1',1,