數(shù)值分析作業(yè)

1、數(shù)值分析實驗作業(yè) 指導老師：院 系： 姓 名： 學 號： 第二章 用matlab編寫的程序如下：function t_charpt2 %數(shù)值實驗二：含“實驗2.1：多項式插值的震蕩現(xiàn)象”和“實驗2.2：樣條插值的收斂”%輸入：實驗選擇，函數(shù)式選擇，插值結(jié)點數(shù)%輸出：擬合函數(shù)及原函數(shù)的圖形result=inputdlg('請選擇實驗，若選2.1，請輸入1，否則輸入2：','charpt2',1,'1');Nb=str2num(char(result);if(Nb=1)&(Nb=2)errordlg('實驗選擇錯誤！');re

2、turn;end promps='請選擇實驗函數(shù)，若選f(x),請輸入f，若選h(x),請輸入h,若選g(x),請輸入g:' titles='charpt2' result=inputdlg(promps,'charpt2',1,'f'); Nb_f=char(result); if(Nb_f='f'&Nb_f='h'&Nb_f='g')errordlg('實驗函數(shù)選擇錯誤！');return;endresult=inputdlg('請輸入插值

3、結(jié)點數(shù)N：','charpt2',1,'10');Nd=str2num(char(result);if(Nd<1)errordlg('結(jié)點輸入錯誤！');return;end switch Nb_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;

4、b=5; end if(Nb=1) 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'); hold on; plot(x,y,'b-'); xlabel('x');ylabel('y=f(x) o and y=Ln(x)-'); elseif(Nb=2) x0=linspace(a,b,Nd+1);y0=feval(f,x0); x=a:0.1:b; cs=spline(x0,y0);y=ppval(cs,x)

5、; plot(x0,y0,'o');hold on ;plot(x,y,'k-'); xlabel('x');ylabel('y=f(x) o and y=Spline(x)-'); end % function y=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; for j=1:n if(j=k) p=p*(z-x0(j)/(x0(k)-x0(j); end end s=s+p*y0(k); end y(

6、i)=s;end若選擇實驗2.1，實驗函數(shù)為f，插值結(jié)點數(shù)為6，則結(jié)果為：若選擇實驗2.1，實驗函數(shù)為f，插值結(jié)點數(shù)為20，則結(jié)果為：若選擇實驗2.1，實驗函數(shù)為f，插值結(jié)點數(shù)為27，則結(jié)果為：選擇其它的函數(shù)重復上述的實驗，在這里我選擇的是h函數(shù)，具體結(jié)果如下：若選擇實驗2.1，實驗函數(shù)為h，插值結(jié)點數(shù)為6，則結(jié)果為：若選擇實驗2.1，實驗函數(shù)為h，插值結(jié)點數(shù)為20，則結(jié)果為：若選擇實驗2.1，實驗函數(shù)為h，插值結(jié)點數(shù)為27，則結(jié)果為：實驗2.2輸入2，實驗函數(shù)為f，插值節(jié)點數(shù)為20，則結(jié)果為：輸入2，實驗函數(shù)為h，插值節(jié)點數(shù)為20，則結(jié)果為：分析：1.對于n次Lagrange插值多項式Pn（

7、x）近似逼近函數(shù)f（x）時，Pn（x）的次數(shù)并非越高，越逼近目標函數(shù)f（x）。在實驗中選取一個較小數(shù)6，可以看到Pn（x）圖像與f（x）吻合情況不是很好，當選取n=20時，Pn（x）與f（x）在區(qū)間-1,1上的吻合情況較好；但是當選取較大數(shù)27，Pn（x）開始發(fā)生振蕩，吻合情況較差，也即“龍格”現(xiàn)象。在選取的函數(shù)h（x）不會發(fā)生類似的現(xiàn)象，精度會隨著次數(shù)的增加而增加，不會出現(xiàn)所謂的“龍格”現(xiàn)象。2.從樣條插值的實驗結(jié)果可以看到：樣條插值的逼近情況比Lagrange插值要好，并且Lagrange插值在區(qū)間-5，5上端點處會有一定的波動，但是樣條插值結(jié)果與精確值吻合較好。第三章 實驗3用matla

8、b編寫的程序如下：function charpt3 %數(shù)值實驗三：含“實驗3.1”和“實驗3.2”%子函數(shù)調(diào)用：dlsa%輸入：實驗選擇%輸出：原函數(shù)及求得的相應(yīng)插值多項式的函數(shù)的圖像以及參數(shù)alph和誤差rresult=inputdlg('請選擇實驗，若選3.1，請輸入1，否則輸入2：','charpt3',1,'1');Nb=str2num(char(result);if(Nb=1)&(Nb=2)errordlg('實驗選擇錯誤！');return;end x0=-1:0.5:2; y0=-4.447 -0.452 0

9、.551 0.048 -0.447 0.549 4.552; if(Nb=1) 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 * and ployfit.y-'); hold on plot(x0,y0,'*') grid on; else result=inputdlg('請

10、輸入權(quán)向量w:','charpt3',1,'1 1 1 1 1 1 1'); w=str2num(char(result); n=3; a,b,c,alph,r=dlsa(x0,y0,w,n); end disp('平方誤差:',num2str(r) disp('參數(shù)alph:',num2str(alph) function a,b,c,alph,r=dlsa(x,y,w,n) %功能：用正交化方法對離散數(shù)據(jù)作多項式最小二乘擬合。 %輸入：m+1個離散點（x,y,w),x,y,w分別用行向量給出。 % 擬合多項式的次數(shù)n,

11、0<n<m. % 平方誤差r=(y-s,y-s),并作離散點列和擬合曲線的圖形 m=length(x)-1; if(n<1|n>=m)errordlg('錯誤：n<1或者n>=m!');return;end%求三項遞推公式的參數(shù)a,b,擬合多項式s(x)的系數(shù)c,其中d(k)=(y,sk); s1=0;s2=ones(1,m+1);v2=sum(w); d(1)=y*w'c(1)=d(1)/v2; for k=1:n xs=x.*s2.2*w'a(k)=xs/v2; if k=1 b(k)=0; else b(k)=v2/v1

12、; end s3=(x-a(k).*s2-b(k)*s1; v3=s3.2*w' d(k+1)=y.*s3*w'c(k+1)=d(k+1)/v3; s1=s2;s2=s3;v1=v2;v2=v3; end %求平方誤差r r=y.*y*w'-c*d'%求擬合多項式s(x)的降冪系數(shù)alph alph=zeros(1,n+1);T=zeros(n+1,n+2); T(:,2)=ones(n+1,1);T(2,3)=-a(1); if n>=2 for k=3:n+1 for i=3:k+1 T(k,i)=T(k-1,i)-a(k-1)*T(k-1,i-1)-

13、b(k-1)*T(k-2,i-2); end end end for i=1:n+1 for k=i:n+1 alph(n+2-i)=alph(n+2-i)+c(k)*T(k,k+2-i); end end %用秦九韶方法計算s(t)的輸出序列（t,s) xmin=min(x);xmax=max(x);dx=(xmax-xmin)/(25*m); t=(xmin-dx):dx:(xmax+dx); s=alph(1); for k=2:n+1 s=s.*t+alph(k); end%輸出點列x-y和擬合曲線t-s的圖形 plot(x,y,'*',t,s,'-')

14、; title('離散數(shù)據(jù)的多項式擬合'); xlabel('x');ylabel('y');grid on;實驗3.1，輸入1，權(quán)函數(shù)為1 1 1 1 1 1 1，結(jié)果為：>> 平方誤差:2.1762e-005參數(shù)alph:1.9991 -2.9977-3.9683e-005 0.54912圖形為：實驗3.2：輸入2，權(quán)函數(shù)為1 1 1 1 1 1 1，結(jié)果為：>> 平方誤差:2.1762e-005參數(shù)alph:1.9991 -2.9977-3.9683e-005 0.54912圖形為：分析：利用最小二乘法作曲線的擬合，

15、對實驗3.1給出的數(shù)據(jù)作的三次多項式的圖形和數(shù)據(jù)節(jié)點擬合得較好。正交化多項式，取w=1 1 1 1 1 1 1時，曲線的擬合也是十分吻合的。由于沒有出現(xiàn)病態(tài)法方程組的情況，二者圖形結(jié)果及平方誤差也相差不大。第四章 實驗4用matlab編寫的程序如下：function t_charpt4%數(shù)值實驗四：含“實驗4.1：復化求積公式計算定積分”和“實驗4.2：高斯數(shù)值積分法用于積分方程求解”%子函數(shù)調(diào)用：CG_L_I function(復化Gauss_Legendre I)公式、CTrapezia(復化梯形公式）、CSimpson(復化Simpson公式）%輸入：實驗選擇、積分法選擇、積分式題號選擇

16、%輸出：積分（實驗4.1）或方程解（實驗4.2）的精確值和數(shù)值解及誤差 result=inputdlg('請選擇實驗,若選4.1,請輸入1,否則輸入2:','charpt4',1,'1'); Nt=str2num(char(result); if(Nt=1)&(Nt=2)errordlg('實驗選擇錯誤！');return;end promps='請選擇積分法,若用復化梯形,輸入T,用復化Simpson,輸入S,用復化Gauss_Legendre,輸入GL:' result=inputdlg(promps,

17、'charpt4',1,'T'); Nb=char(result); if(Nb='T' & Nb='S' & Nb='GL')errordlg('積分公式選擇錯誤！');return;end if(Nt=1) result=inputdlg('請輸入積分式題號14:','實驗4.1',1,'1'); Nb_f=str2num(char(result); if(Nb_f<1)|(Nb_f>4)errordlg('沒有

18、該積分式!');return;end switch Nb_f case 1 fun=inline('-2./(x.2-1)');a=2;b=3; case 2 fun=inline('4./(x.2+1)');a=0;b=1; case 3 fun=inline('3.x');a=0;b=1; case 4 fun=inline('x.*exp(x)');a=1;b=2; end tol=0.5e-7;h=0.01; if(Nb='T')%用復化梯形公式 t=(fun(a)+fun(b)*(b-a)/2;

19、k=1;t0=0; while(abs(t-t0)>=tol*3) t0=t;h=(b-a)/2k; t=t0/2+h*sum(fun(a+h:2*h:b-h); k=k+1; end elseif(Nb='S')%用復化Simpson公式 t=quad(fun,a,b,tol); elseif(Nb='GL')%用復化Gauss_Legendre I N=floor(b-a)/h);t=0;xk=0; for k=0:N xk=a+k*h+h/2; t=t+fun(xk-h/(2*sqrt(3)+fun(xk+h/(2*sqrt(3); end t=t*

20、h/2; end elseif(Nt=2) result=inputdlg('請輸入方程式題號1或2:','實驗4.2',1,'1'); Nb_f=str2num(char(result); if(Nb_f=1 & Nb_f=2)errordlg('沒有該方程式!');return;end result=inputdlg('請輸入步長:','實驗4.2',1,'0.01'); h=str2num(char(result); if(h<=0)errordlg('請

21、輸入正確的步長!');return;end if(Nb='T')%用復化梯形公式 x,t=CTrapezia(0,1,h,Nb_f); elseif(Nb='S')%用復化Simpson公式 x,t=CSimpson(0,1,h,Nb_f); elseif(Nb='GL')%用復化Gauss_Legendre I公式 x,t=CG_L_I(0,1,h,Nb_f); end plot(x,t,'g-'); xlabel('x');ylabel('y'); title('積分方程求解&#

22、39;) hold on disp('實驗4.2(',num2str(Nb_f),')的計算結(jié)果:',num2str(t'); if(Nb_f=1) fplot('exp(x)',0 1,'*'); hold off disp('實驗4.2題（1）的精確解:',num2str(exp(x'); disp('絕對誤差和:',num2str(sum(abs(t'-exp(x'); else fplot('1/(1+t)2',0 1,'*')

23、; hold off y=1./(x+1).2; disp('實驗4.2題（2）的精確解:',num2str(exp(y'); disp('絕對誤差和:',num2str(sum(abs(t'-y'); end end if(Nt=1) disp('實驗4.1題(',num2str(Nb_f),')的計算結(jié)果:',num2str(t); switch Nb_f case 1 disp('精確解:ln2-ln3=-0.4054651081') disp('絕對誤差:',num2

24、str(abs(t+0.4054651081); case 2 disp('精確解:pi=3.14159265358979') disp('絕對誤差:',num2str(abs(t-pi); case 3 disp('精確解:2/ln3=1.82047845325368') disp('絕對誤差:',num2str(abs(t-1.82047845325368); case 4 disp('精確解:e2=7.38905609893065') disp('絕對誤差:',num2str(abs(t-7.

25、38905609893065); endend %function x,y=CG_L_I(a,b,h,N)%復化Gauss_Legendre I,用于積分方程求解%輸入：a、b分別為求積下、上限，h為步長，N為方程式題號實驗選擇、積分法選擇、積分式題號選擇%輸出：x,y為方程離散解y(i)=f(x(i) n=floor(b-a)/h); A=zeros(2*n+2,2*n+2);b=zeros(2*n+2,1);x=b; for i=0:n t1=a+h*(i+0.5-sqrt(3)/6);t2=a+h*(i+0.5+sqrt(3)/6); x(2*i+1)=t1;x(2*i+2)=t2; i

26、f(N=1) b(2*i+1)=exp(t1);b(2*i+2)=exp(t2); else b(2*i+1)=-(4*t1.3+5*t1.2-2*t1+5)/(8*(t1+1).2); b(2*i+2)=-(4*t2.3+5*t2.2-2*t2+5)/(8*(t2+1).2); end for j=0:n x1=a+h*(j+0.5-sqrt(3)/6);x2=a+h*(j+0.5+sqrt(3)/6); if(N=1) A(2*i+1,2*j+1)=exp(t1)*2/(exp(1)-1); A(2*i+1,2*j+2)=exp(t1)*2/(exp(1)-1); A(2*i+2,2*j+

27、1)=exp(t2)*2/(exp(1)-1); A(2*i+2,2*j+2)=exp(t2)*2/(exp(1)-1); else A(2*i+1,2*j+1)=1/(1+x1)-t1; A(2*i+1,2*j+2)=1/(1+x2)-t1; A(2*i+2,2*j+1)=1/(1+x1)-t2; A(2*i+2,2*j+2)=1/(1+x2)-t2; end end endA=h*A/2-eye(2*n+2);y=inv(A)*b; %function x,y=CTrapezia(a,b,h,N)%復化Gauss_%復化梯形公式，用于積分方程求解%輸入：a、b分別為求積下、上限，h為步長，

28、N為方程式題號實驗選擇、積分法選擇、積分式題號選擇%輸出：x,y為方程離散解y(i)=f(x(i) n=floor(b-a)/h); A=zeros(n+1,n+1);b=zeros(n+1,1);x=b; for i=0:n t=a+i*h;x(i+1)=t; if(N=1) b(i+1)=-exp(t); else b(i+1)=(4*t.3+5*t.2-2*t+5)/(8*(t+1).2); end for j=0:n s=a+j*h; if(j=0)|(j=n) if(N=1) A(i+1,j+1)=exp(t)*2/(exp(1)-1); else A(i+1,j+1)=1/(1+s

29、)-t; end else if(N=1) A(i+1,j+1)=exp(t)*2/(exp(1)-1)*2; else A(i+1,j+1)=(1/(1+s)-t)*2; end end end endA=A-eye(n+1)*2/h;b=-2*b/h;y=inv(A)*b; %function x,y=CSimpson(a,b,h,N)%復化Simpson公式，用于積分方程求解%輸入：a、b分別為求積下、上限，h為步長，N為方程式題號實驗選擇、積分法選擇、積分式題號選擇%輸出：x,y為方程離散解y(i)=f(x(i) n=floor(b-a)/h);h1=h/2; A=zeros(2*n+

30、1,2*n+1);b=zeros(2*n+1,1);x=b; for i=0:n*2 t=a+i*h1;x(i+1)=t; if(N=1) b(i+1)=-exp(t); else b(i+1)=(4*t.3+5*t.2-2*t+5)/(8*(t+1).2); end for j=0:n*2 s=a+j*h1; if(j=0)|(j=n*2) if(N=1) A(i+1,j+1)=exp(t)*2/(exp(1)-1); else A(i+1,j+1)=1/(1+s)-t; end elseif(mod(j,2)=0) if(N=1) A(i+1,j+1)=exp(t)*2/(exp(1)-1

31、)*2; else A(i+1,j+1)=(1/(1+s)-t)*2; end elseif(mod(j,2)=0) if(N=1) A(i+1,j+1)=exp(t)*2/(exp(1)-1)*4; else A(i+1,j+1)=(1/(1+s)-t)*4; end end end end A=A-eye(2*n+1)*6/h;b=-6*b/h; y=inv(A)*b;1.若選擇實驗4.1，1）積分法選擇T(即復化梯形公式)：>> 實驗4.1題(1)的計算結(jié)果:-0.40547精確解:ln2-ln3=-0.4054651081絕對誤差:1.3944e-008>> 實

32、驗4.1題(2)的計算結(jié)果:3.1416精確解:pi=3.14159265358979絕對誤差:3.9736e-008>> 實驗4.1題(3)的計算結(jié)果:1.8205精確解:2/ln3=1.82047845325368絕對誤差:4.3655e-008>> 實驗4.1題(4)的計算結(jié)果:7.3891精確解:e2=7.38905609893065絕對誤差:2.0775e-0082）積分法選擇S(即復化Simpson公式)：>> 實驗4.1題(1)的計算結(jié)果:-0.40547精確解:ln2-ln3=-0.4054651081絕對誤差:1.2625e-009>

33、> 實驗4.1題(2)的計算結(jié)果:3.1416精確解:pi=3.14159265358979絕對誤差:2.7517e-010>> 實驗4.1題(3)的計算結(jié)果:1.8205精確解:2/ln3=1.82047845325368絕對誤差:1.0877e-010>> 實驗4.1題(4)的計算結(jié)果:7.3891精確解:e2=7.38905609893065絕對誤差:8.2168e-0113）積分法選擇GL(即復化Gauss_Legendre公式):>> 實驗4.1題(1)的計算結(jié)果:-0.40796精確解:ln2-ln3=-0.4054651081絕對誤差:0

34、.0024907>> 實驗4.1題(2)的計算結(jié)果:3.1615精確解:pi=3.14159265358979絕對誤差:0.0199>> 實驗4.1題(3)的計算結(jié)果:1.8506精確解:2/ln3=1.82047845325368絕對誤差:0.030165>> 實驗4.1題(4)的計算結(jié)果:7.538精確解:e2=7.38905609893065絕對誤差:0.148892.若選擇實驗4.2，1）積分法選擇T(即復化梯形公式),題號選擇1，積分步長為0.01，則結(jié)果如下：>> 實驗4.2(1)的計算結(jié)果:0.99998 1.01 1.0202 1

35、.0304 1.0408 1.0513 1.0618 1.0725 1.0833 1.0942 1.1052 1.1163 1.1275 1.1388 1.1503 1.1618 1.1735 1.1853 1.1972 1.2092 1.2214 1.2337 1.2461 1.2586 1.2712 1.284 1.2969 1.3099 1.3231 1.3364 1.3498 1.3634 1.3771 1.3909 1.4049 1.419 1.4333 1.4477 1.4623 1.477 1.4918 1.5068 1.5219 1.5372 1.5527 1.5683 1.5

36、84 1.6 1.616 1.6323 1.6487 1.6653 1.682 1.6989 1.716 1.7332 1.7506 1.7682 1.786 1.804 1.8221 1.8404 1.8589 1.8776 1.8964 1.9155 1.9348 1.9542 1.9738 1.9937 2.0137 2.034 2.0544 2.075 2.0959 2.117 2.1382 2.1597 2.1814 2.2034 2.2255 2.2479 2.2705 2.2933 2.3163 2.3396 2.3631 2.3869 2.4109 2.4351 2.4596

37、2.4843 2.5092 2.5345 2.5599 2.5857 2.6117 2.6379 2.6644 2.6912 2.7182實驗4.2題（1）的精確解:1 1.0101 1.0202 1.0305 1.0408 1.0513 1.0618 1.0725 1.0833 1.0942 1.1052 1.1163 1.1275 1.1388 1.1503 1.1618 1.1735 1.1853 1.1972 1.2092 1.2214 1.2337 1.2461 1.2586 1.2712 1.284 1.2969 1.31 1.3231 1.3364 1.3499 1.3634 1

38、.3771 1.391 1.4049 1.4191 1.4333 1.4477 1.4623 1.477 1.4918 1.5068 1.522 1.5373 1.5527 1.5683 1.5841 1.6 1.6161 1.6323 1.6487 1.6653 1.682 1.6989 1.716 1.7333 1.7507 1.7683 1.786 1.804 1.8221 1.8404 1.8589 1.8776 1.8965 1.9155 1.9348 1.9542 1.9739 1.9937 2.0138 2.034 2.0544 2.0751 2.0959 2.117 2.138

39、3 2.1598 2.1815 2.2034 2.2255 2.2479 2.2705 2.2933 2.3164 2.3396 2.3632 2.3869 2.4109 2.4351 2.4596 2.4843 2.5093 2.5345 2.56 2.5857 2.6117 2.6379 2.6645 2.6912 2.7183絕對誤差和:0.00289482）積分法選擇S(即復化Simpson公式),題號選擇1，積分步長為0.01，則結(jié)果如下>> 實驗4.2(1)的計算結(jié)果:1 1.005 1.0101 1.0151 1.0202 1.0253 1.0305 1.0356 1

40、.0408 1.046 1.0513 1.0565 1.0618 1.0672 1.0725 1.0779 1.0833 1.0887 1.0942 1.0997 1.1052 1.1107 1.1163 1.1219 1.1275 1.1331 1.1388 1.1445 1.1503 1.156 1.1618 1.1677 1.1735 1.1794 1.1853 1.1912 1.1972 1.2032 1.2092 1.2153 1.2214 1.2275 1.2337 1.2399 1.2461 1.2523 1.2586 1.2649 1.2712 1.2776 1.284 1.2

41、905 1.2969 1.3034 1.31 1.3165 1.3231 1.3298 1.3364 1.3431 1.3499 1.3566 1.3634 1.3703 1.3771 1.384 1.391 1.3979 1.4049 1.412 1.4191 1.4262 1.4333 1.4405 1.4477 1.455 1.4623 1.4696 1.477 1.4844 1.4918 1.4993 1.5068 1.5144 1.522 1.5296 1.5373 1.545 1.5527 1.5605 1.5683 1.5762 1.5841 1.592 1.6 1.608 1.

42、6161 1.6242 1.6323 1.6405 1.6487 1.657 1.6653 1.6736 1.682 1.6905 1.6989 1.7074 1.716 1.7246 1.7333 1.7419 1.7507 1.7594 1.7683 1.7771 1.786 1.795 1.804 1.813 1.8221 1.8313 1.8404 1.8497 1.8589 1.8682 1.8776 1.887 1.8965 1.906 1.9155 1.9251 1.9348 1.9445 1.9542 1.964 1.9739 1.9838 1.9937 2.0037 2.0138 2.0238 2.034 2.0442 2.0544 2.0647 2.0751 2.0855 2.0959 2.1064 2.117 2.1276 2.1383 2.149 2.1598 2.1706 2.1815 2.1924 2.2034 2.2144 2.2255 2.2367 2.2479 2.2592 2.2705 2.2819 2.2933 2.3048 2.3164 2.328 2.3396 2.3514 2.3632 2.375 2.3869 2.3989 2.4109 2.423