數(shù)值分析上機(jī)題目

1、實(shí)驗(yàn)一實(shí)驗(yàn)項(xiàng)目：共軛梯度法求解對(duì)稱正定的線性方程組實(shí)驗(yàn)內(nèi)容：用共軛梯度法求解下面方程組(1) 迭代20次或滿足時(shí)停止計(jì)算。編制程序：儲(chǔ)存m文件function x,k=CGmethod(A,b)n=length(A);x=2*ones(n,1);r=b-A*x;rho=r*r;k=0;while rho10(-11) & k10(-7)&k104 k=k+1; if k=1 p=r; else beta=(r1*r1)/(r*r);p=r1+beta*p; end r=r1; w=A*p; alpha=(r*r)/(p*w); x=x+alpha*p; r1=r-alpha*w;end運(yùn)行程序

2、：clear,clcn=1000;A=hilb(n);b=sum(A);x,k=CGmethod(A,b)實(shí)驗(yàn)二1、 實(shí)驗(yàn)?zāi)康模河脧?fù)化Simpson方法、自適應(yīng)復(fù)化梯形方法和Romberg方法求數(shù)值積分。實(shí)驗(yàn)內(nèi)容：計(jì)算下列定積分(1) (2) (3) 實(shí)驗(yàn)要求：（1）分別用復(fù)化Simpson公式、自適應(yīng)復(fù)化梯形公式計(jì)算要求絕對(duì)誤差限為，輸出每種方法所需的節(jié)點(diǎn)數(shù)和積分近似值，對(duì)于自適應(yīng)方法，顯示實(shí)際計(jì)算節(jié)點(diǎn)上離散函數(shù)值的分布圖；（2）分析比較計(jì)算結(jié)果。程序：syms xf=x6/10-x2+x %定義函數(shù)f（x）n=input(輸入所求導(dǎo)數(shù)階數(shù):) f2=diff(f,x,n) %求f(x)的n

3、階導(dǎo)數(shù)（1）復(fù)化梯形clcclearsyms x %定義自變量xf=inline(x6/10-x2+x,x) %定義函數(shù)f(x)=x*exp(x)，換函數(shù)時(shí)只需換該函數(shù)表達(dá)式即可f2=inline(3*x4 - 2,x) %定義f(x)的二階導(dǎo)數(shù),輸入程序1里求出的f2即可。f3=-(3*x4 - 2) %因fminbnd（）函數(shù)求的是表達(dá)式的最小值，且要求表達(dá)式帶引號(hào)，故取負(fù)號(hào)，以便求最大值e=0.5*10(-7) %精度要求值 a=0 %積分下限b=2 %積分上限x1=fminbnd(f3,1,2) %求負(fù)的二階導(dǎo)數(shù)的最小值點(diǎn)，也就是求二階導(dǎo)數(shù)的最大值點(diǎn)對(duì)應(yīng)的x值for n=2: %求等分

4、數(shù)n Rn=-(b-a)/12*(b-a)/n)2*f2(x1) %計(jì)算余項(xiàng) if abs(Rn)e %用余項(xiàng)進(jìn)行判斷 break % 符合要求時(shí)結(jié)束 endendh=(b-a)/n %求hTn1=0 for k=1:n-1 %求連加和 xk=a+k*h Tn1=Tn1+f(xk)endTn=h/2*(f(a)+2*Tn1+f(b)z=exp(2)R=Tn-z %求已知值與計(jì)算值的差stem(xk,Tn1);fprintf(用復(fù)化梯形算法計(jì)算的結(jié)果 Tn=)disp(Tn)fprintf(等分?jǐn)?shù) n=)disp(n) %輸出等分?jǐn)?shù)fprintf(已知值與計(jì)算值的誤差 R=)disp(R)（2）

5、復(fù)化Simpsonclcclearsyms x %定義自變量xf=inline(x6/10-x2+x,x) %定義函數(shù)f(x)=x*exp(x)，換函數(shù)時(shí)只需換該函數(shù)表達(dá)式即可f2=inline(36*x2,x) %定義f(x)的四階導(dǎo)數(shù)，輸入程序1里求出的f2即可f3=-(36*x2) %因fminbnd（）函數(shù)求的是表達(dá)式的最小值，且要求表達(dá)式帶引號(hào)，故取負(fù)號(hào)，一邊求最大值e=5*10(-8) %精度要求值 a=0 %積分下限b=2 %積分上限x1=fminbnd(f3,1,2) %求負(fù)的四階導(dǎo)數(shù)的最小值點(diǎn)，也就是求四階導(dǎo)數(shù)的最大值點(diǎn)對(duì)應(yīng)的x值for n=2: %求等分?jǐn)?shù)n Rn=-(b-

6、a)/180*(b-a)/(2*n)4*f2(x1) %計(jì)算余項(xiàng) if abs(Rn)e %用余項(xiàng)進(jìn)行判斷 break % 符合要求時(shí)結(jié)束 endendh=(b-a)/n %求hSn1=0 Sn2=0for k=0:n-1 %求兩組連加和 xk=a+k*h xk1=xk+h/2 Sn1=Sn1+f(xk1) Sn2=Sn2+f(xk)end Sn=h/6*(f(a)+4*Sn1+2*(Sn2-f(a)+f(b) %因Sn2多加了k=0時(shí)的值，故減去f(a)z=exp(2)R=Sn-z %求已知值與計(jì)算值的差fprintf(用Simpson公式計(jì)算的結(jié)果 Sn=)disp(Sn)fprintf(

7、等分?jǐn)?shù) n=)disp(n) fprintf(已知值與計(jì)算值的誤差 R=)disp(R)結(jié)果：用復(fù)化梯形算法計(jì)算的結(jié)果 Tn= 1.6674等分?jǐn)?shù) n= 24764已知值與計(jì)算值的誤差 R= -6.3976用Simpson公式計(jì)算的結(jié)果 Sn= 1.0434等分?jǐn)?shù) n= 76已知值與計(jì)算值的誤差 R= -6.0216用復(fù)化梯形算法計(jì)算的結(jié)果 Tn= 0.8985等分?jǐn)?shù) n= 1119已知值與計(jì)算值的誤差 R= -6.1665用Simpson公式計(jì)算的結(jié)果 Sn= 0.9406等分?jǐn)?shù) n= 8已知值與計(jì)算值的誤差 R= -6.1245用復(fù)化梯形算法計(jì)算的結(jié)果 Tn= 23.2353等分?jǐn)?shù) n=

8、已知值與計(jì)算值的誤差 R= 16.1704用Simpson公式計(jì)算的結(jié)果 Sn= 23.2617等分?jǐn)?shù) n= 10647已知值與計(jì)算值的誤差 R= 16.1969分析：在處理問(wèn)題時(shí)，復(fù)化Simpson要比復(fù)化梯度計(jì)算速度要快很多。2、實(shí)驗(yàn)?zāi)康模焊咚箶?shù)值積分方法用于積分方程求解。實(shí)驗(yàn)內(nèi)容：線性的積分方程的數(shù)值求解，可以被轉(zhuǎn)化為線性代數(shù)方程組的求解問(wèn)題。而線性代數(shù)方程組所含未知數(shù)的個(gè)數(shù)，與用來(lái)離散積分的數(shù)值方法的節(jié)點(diǎn)個(gè)數(shù)相同。在節(jié)點(diǎn)數(shù)相同的前提下，高斯數(shù)值積分方法有較高的代數(shù)精度，用它通常會(huì)得到較好的結(jié)果。對(duì)第二類(lèi)Fredholm積分方程首先將積分區(qū)間a,b等分成n份，在每個(gè)子區(qū)間上離散方程中的積

9、分就得到線性代數(shù)方程組。實(shí)驗(yàn)要求：分別使用如下方法，離散積分方程中的積分1.復(fù)化梯形方法；2.復(fù)化辛甫森方法；3.復(fù)化高斯方法。求解如下的積分方程，方程的準(zhǔn)確解為，并比較各算法的優(yōu)劣。程序結(jié)果：當(dāng)?shù)螖?shù)n=1時(shí)精確解1.0000 1.2840 1.6487 2.1170 2.7183 復(fù)化梯形方法0.8591 1.1032 1.4165 1.81882.3354復(fù)化辛甫森方法0.9993 1.2832 1.6476 2.1156 2.7165復(fù)化高斯方法1.0004 1.2846 1.6495 2.1180 2.7195復(fù)化梯形方法的平均誤差err=0.247復(fù)化辛甫森方法的平均誤差err=

10、0.00116復(fù)化高斯方法的平均誤差err=0.0008當(dāng)?shù)螖?shù)n=5時(shí)，精確解1.0000 1.2840 1.6487 2.1170 2.7183 復(fù)化梯形方法0.9934 1.2755 1.6378 2.1030 2.7003 復(fù)化辛甫森方法1.0000 1.2840 1.6487 2.1170 2.7183復(fù)化高斯方法1.0000 1.2840 1.6487 2.1170 2.7183復(fù)化梯形方法的平均誤差err=0.00116復(fù)化辛甫森方法和復(fù)化高斯方法的平均誤差err=0可以看出，復(fù)化高斯方法得到的結(jié)果精度最高，復(fù)化辛普森方法比復(fù)化梯形方法的精度要高，程序：clear,clca=0

11、;b=1;n=5;figurefun1(a,b,n);hold ona=0;b=1;n=5;figurefun2(a,b,n);hold onfigurefun3(a,b,n);編制m函數(shù)：function y=Kfun(t,s)y=2/(exp(1)-1)*exp(t);function y=ffun(t)y=-exp(t);function y=Fexc(t)%精確解y=exp(t);function x,w=fhgauss(a,b,n)h=(b-a)/n;x1=a:h:b;x=zeros(1,2*n);w=x;for i=1:nx(2*i-1:2*i),w(2*i-1:2*i)=Gaus

12、sLegendre(x1(i),x1(i+1),2);endfunction x,w=fhsimpson(a,b,n)h=(b-a)/n;x=a:h/2:b;w=x;w(1)=h/6;w(2*n+1)=h/6;for i=1:nw(2*i)=2/3*h;if in w(2*i+1)=h/3;endendfunction x,w=fhtrapz(a,b,n)h=(b-a)/n;x=a:h:b;for i=1:n+1 if i=1|i=n+1 w(i)=h/2; else w(i)=h; endendfunction x,w=GaussLegendre(a,b,n)i=1:n-1;c=i./sqr

13、t(4*i.2-1);CM=diag(c,1) + diag(c,-1);V L=eig(CM);x ind=sort(diag(L);V=V(:,ind);w=2*V(:,1).2;x=x*(b-a)/2+(b+a)/2;w=w*(b-a)/2;function fun1(a,b,n)x1,w1=fhtrapz(a,b,n);n1=4;n=n+1;h=(b-a)/n1;x=a:h:b;y0=Fexc(x);A=zeros(n,n);B=zeros(n,1);for i=1:nB(i)=ffun(x1(i);endfor i=1:nfor j=1:n A(i,j)=w1(j)*Kfun(x1(

14、i),x1(j); endendA=eye(n)-A;y1=(AB);yN=x;for i=1:n1+1yN(i)=ffun(x(i);for j=1:nyN(i)=yN(i)+w1(j)*Kfun(x(i),x1(j)*y1(j);endendfprintf(數(shù)值解:n)yNfprintf(精確解:n)y0plot(x,yN,x,x,y0,.);h=legend(復(fù)化值,真實(shí)值);function fun2(a,b,n)x1,w1=fhsimpson(a,b,n);n=2*n+1;n1=50;h=(b-a)/n1;x=a:h:b;y0=Fexc(x);A=zeros(n,n);B=zeros

15、(n,1);for i=1:nB(i)=ffun(x1(i);endfor i=1:nfor j=1:n A(i,j)=w1(j)*Kfun(x1(i),x1(j); endendA=eye(n)-A;y1=(AB);yN=x;for i=1:n1+1yN(i)=ffun(x(i);for j=1:nyN(i)=yN(i)+w1(j)*Kfun(x(i),x1(j)*y1(j);endendfprintf(數(shù)值解:n)yNfprintf(精確解:n)y0plot(x,yN,x,x,y0,.);h=legend(復(fù)化值,真實(shí)值);function fun3(a,b,n)x1,w1=fhgauss

16、(a,b,n);n=2*n;n1=4;h=(b-a)/n1;x=a:h:b;y0=Fexc(x);A=zeros(n,n);B=zeros(n,1);for i=1:nB(i)=ffun(x1(i);endfor i=1:nfor j=1:n A(i,j)=w1(j)*Kfun(x1(i),x1(j); endendA=eye(n)-A;y1=(AB);yN=x;for i=1:n1+1yN(i)=ffun(x(i);for j=1:nyN(i)=yN(i)+w1(j)*Kfun(x(i),x1(j)*y1(j);endendfprintf(數(shù)值解:n)yNfprintf(精確解:n)y0pl

17、ot(x,yN,x,x,y0,.);h=legend(復(fù)化值,真實(shí)值);圖一圖二圖三實(shí)驗(yàn)三1、對(duì)常微分方程初值問(wèn)題 分別使用Euler顯示方法、改進(jìn)的Euler方法和經(jīng)典RK法和四階Adams預(yù)測(cè)-校正算法，求解常微分方程數(shù)值解，并與其精確解進(jìn)行作圖比較。其中多步法需要的初值由經(jīng)典RK法提供。程序：子程序function E=Euler(fun,x0,y0,h,N)x=zeros(1,N+1);y=zeros(1,N+1);x(1)=x0;y(1)=y0;for n=1:Nx(n+1)=x(n)+h;y(n+1)=y(n)+h*feval(fun,x(n),y(n);endx1=xy1=ypl

18、ot(x1,y1,k)hold onfunction E=Euler_modify(fun,x0,y0,h,N)x=zeros(1,N+1);y=zeros(1,N+1);x(1)=x0;y(1)=y0;for n=1:Nx(n+1)=x(n)+h;z0=y(n)+h*feval(fun,x(n),y(n);y(n+1)=y(n)+h/2*(feval(fun,x(n),y(n)+feval(fun,x(n+1),z0);endx2=xy2=yplot(x2,y2,g)function x,y=Rk_N4(f,x0,y0,h,N) x=zeros(1,N+1);y=zeros(1,N+1);

19、x(1)=x0;y(1)=y0; for n=1:N x(n+1)=x(n)+h; k1=h*feval(f,x(n),y(n); k2=h*feval(f,x(n)+1/2*h,y(n)+1/2*k1); k3=h*feval(f,x(n)+1/2*h,y(n)+1/2*k2); k4=h*feval(f,x(n)+h,y(n)+k3); y(n+1)=y(n)+1/6*(k1+2*k2+2*k3+k4); end 運(yùn)行以下程序clcclearEuler(fun,0,1,0.1,100)Euler_modify(fun,0,1,0.1,100)x,y=Rk_N4(fun,0,1,0.1,10

20、0)plot(x,y,-*r)hold onx1=0:0.1:pi;y1=cos(x1)+sin(x1)plot(x1,y1,-.b)title(誤差分析);xlabel(x軸);ylabel(y軸);legend(Euler,Euler改進(jìn)法,R_K法,精確);axis(0 pi -1.5 1.5);grid on畫(huà)出圖形進(jìn)行比較：2、實(shí)驗(yàn)?zāi)康模篖orenz問(wèn)題與混沌實(shí)驗(yàn)內(nèi)容：考慮著名的Lorenz方程 其中s, r, b為變化區(qū)域有一定限制的實(shí)參數(shù)。該方程形式簡(jiǎn)單，表面上看并無(wú)驚人之處，但由該方程揭示出的許多現(xiàn)象，促使“混沌”成為數(shù)學(xué)研究的嶄新領(lǐng)域，在實(shí)際應(yīng)用中也產(chǎn)生了巨大的影響。實(shí)驗(yàn)方法：先取定初值Y0=(x, y, z)=(0, 0, 0)，參數(shù)s=10, r=28, b=8/3，用MATLAB編程數(shù)值求解，并與MATLAB函數(shù)ods45的計(jì)算結(jié)果進(jìn)行對(duì)比。實(shí)驗(yàn)要求：（1）對(duì)目前取定的參數(shù)值s, r和b，選取不同的初值Y0進(jìn)行運(yùn)算，觀察計(jì)算的結(jié)果有什么特點(diǎn)？解的曲線是否有界？解的曲線是不是周期的或趨于某個(gè)固定點(diǎn)？（2）在問(wèn)題允許的范圍內(nèi)適當(dāng)改變其中的參數(shù)值s, r, b，再選取不同的初始值Y0進(jìn)行試算，觀察并記錄計(jì)算的結(jié)果有什么特點(diǎn)？是否發(fā)現(xiàn)什么不同的現(xiàn)象？