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

--本頁僅作為文檔封面，使用時請直接刪除即可--

--內(nèi)頁可以根據(jù)需求調(diào)整合適字體及大小本頁僅作為文檔封面，使用時請直接刪除即可--

--內(nèi)頁可以根據(jù)需求調(diào)整合適字體及大小--數(shù)值分析上機(jī)題目4(總21頁)PAGE實(shí)驗(yàn)一實(shí)驗(yàn)項目：共軛梯度法求解對稱正定的線性方程組實(shí)驗(yàn)內(nèi)容：用共軛梯度法求解下面方程組(1)迭代20次或滿足時停止計算。編制程序：儲存m文件function[x,k]=CGmethod(A,b)n=length(A);x=2*ones(n,1);r=b-A*x;rho=r'*r;k=0;whilerho>10^(-11)&k<1000k=k+1;ifk==1p=r;elsebeta=rho/rho1;p=r+beta*p;endw=A*p;alpha=rho/(p'*w);x=x+alpha*p;r=r-alpha*w;rho1=rho;rho=r'*r;end運(yùn)行程序：clear,clcA=[2-100;-13-10;0-14-1;00-15];b=[3-215]';[x,k]=CGmethod(A,b)運(yùn)行結(jié)果：x=(2),A是1000階的Hilbert矩陣或如下的三對角矩陣，A[i,i]=4，A[i,i-1]=A[i-1,i]=-1，i=2,3,..,nb[1]=3,b[n]=3,b[i]=2，i=2,3,…,n-1迭代10000次或滿足時停止計算。編制程序：儲存m文件function[x,k]=CGmethod_1(A,b)n=length(A);x(1:n,1)=0;r=b-A*x;r1=r;k=0;whilenorm(r1,1)>10^(-7)&k<10^4k=k+1;ifk==1p=r;elsebeta=(r1'*r1)/(r'*r);p=r1+beta*p;endr=r1;w=A*p;alpha=(r'*r)/(p'*w);x=x+alpha*p;r1=r-alpha*w;end運(yùn)行程序：clear,clcn=1000;A=hilb(n);b=sum(A')';[x,k]=CGmethod(A,b)實(shí)驗(yàn)二實(shí)驗(yàn)?zāi)康模河脧?fù)化Simpson方法、自適應(yīng)復(fù)化梯形方法和Romberg方法求數(shù)值積分。實(shí)驗(yàn)內(nèi)容：計算下列定積分(1)(2)(3)實(shí)驗(yàn)要求：（1）分別用復(fù)化Simpson公式、自適應(yīng)復(fù)化梯形公式計算要求絕對誤差限為，輸出每種方法所需的節(jié)點(diǎn)數(shù)和積分近似值，對于自適應(yīng)方法，顯示實(shí)際計算節(jié)點(diǎn)上離散函數(shù)值的分布圖；（2）分析比較計算結(jié)果。程序：symsxf=x^6/10-x^2+x%定義函數(shù)f（x）n=input('輸入所求導(dǎo)數(shù)階數(shù):')f2=diff(f,x,n)%求f(x)的n階導(dǎo)數(shù)（1）復(fù)化梯形clcclearsymsx%定義自變量xf=inline('x^6/10-x^2+x','x')%定義函數(shù)f(x)=x*exp(x)，換函數(shù)時只需換該函數(shù)表達(dá)式即可f2=inline('3*x^4-2','x')%定義f(x)的二階導(dǎo)數(shù),輸入程序1里求出的f2即可。f3='-(3*x^4-2)'%因fminbnd（）函數(shù)求的是表達(dá)式的最小值，且要求表達(dá)式帶引號，故取負(fù)號，以便求最大值e=*10^(-7)%精度要求值a=0%積分下限b=2%積分上限x1=fminbnd(f3,1,2)%求負(fù)的二階導(dǎo)數(shù)的最小值點(diǎn)，也就是求二階導(dǎo)數(shù)的最大值點(diǎn)對應(yīng)的x值forn=2:1000000%求等分?jǐn)?shù)nRn=-(b-a)/12*((b-a)/n)^2*f2(x1)%計算余項ifabs(Rn)<e%用余項進(jìn)行判斷break%符合要求時結(jié)束endendh=(b-a)/n%求hTn1=0fork=1:n-1%求連加和xk=a+k*hTn1=Tn1+f(xk)endTn=h/2*((f(a)+2*Tn1+f(b)))z=exp(2)R=Tn-z%求已知值與計算值的差stem(xk,Tn1);fprintf('用復(fù)化梯形算法計算的結(jié)果Tn=')disp(Tn)fprintf('等分?jǐn)?shù)n=')disp(n)%輸出等分?jǐn)?shù)fprintf('已知值與計算值的誤差R=')disp(R)（2）復(fù)化Simpsonclcclearsymsx%定義自變量xf=inline('x^6/10-x^2+x','x')%定義函數(shù)f(x)=x*exp(x)，換函數(shù)時只需換該函數(shù)表達(dá)式即可f2=inline('36*x^2','x')%定義f(x)的四階導(dǎo)數(shù)，輸入程序1里求出的f2即可f3='-(36*x^2)'%因fminbnd（）函數(shù)求的是表達(dá)式的最小值，且要求表達(dá)式帶引號，故取負(fù)號，一邊求最大值e=5*10^(-8)%精度要求值a=0%積分下限b=2%積分上限x1=fminbnd(f3,1,2)%求負(fù)的四階導(dǎo)數(shù)的最小值點(diǎn)，也就是求四階導(dǎo)數(shù)的最大值點(diǎn)對應(yīng)的x值forn=2:1000000%求等分?jǐn)?shù)nRn=-(b-a)/180*((b-a)/(2*n))^4*f2(x1)%計算余項ifabs(Rn)<e%用余項進(jìn)行判斷break%符合要求時結(jié)束endendh=(b-a)/n%求hSn1=0Sn2=0fork=0:n-1%求兩組連加和xk=a+k*hxk1=xk+h/2Sn1=Sn1+f(xk1)Sn2=Sn2+f(xk)endSn=h/6*(f(a)+4*Sn1+2*(Sn2-f(a))+f(b))%因Sn2多加了k=0時的值，故減去f(a)z=exp(2)R=Sn-z%求已知值與計算值的差fprintf('用Simpson公式計算的結(jié)果Sn=')disp(Sn)fprintf('等分?jǐn)?shù)n=')disp(n)fprintf('已知值與計算值的誤差R=')disp(R)結(jié)果：用復(fù)化梯形算法計算的結(jié)果Tn=等分?jǐn)?shù)n=24764已知值與計算值的誤差R=用Simpson公式計算的結(jié)果Sn=等分?jǐn)?shù)n=76已知值與計算值的誤差R=用復(fù)化梯形算法計算的結(jié)果Tn=等分?jǐn)?shù)n=1119已知值與計算值的誤差R=用Simpson公式計算的結(jié)果Sn=等分?jǐn)?shù)n=8已知值與計算值的誤差R=用復(fù)化梯形算法計算的結(jié)果Tn=等分?jǐn)?shù)n=1000000已知值與計算值的誤差R=用Simpson公式計算的結(jié)果Sn=等分?jǐn)?shù)n=10647已知值與計算值的誤差R=分析：在處理問題時，復(fù)化Simpson要比復(fù)化梯度計算速度要快很多。2、實(shí)驗(yàn)?zāi)康模焊咚箶?shù)值積分方法用于積分方程求解。實(shí)驗(yàn)內(nèi)容：線性的積分方程的數(shù)值求解，可以被轉(zhuǎn)化為線性代數(shù)方程組的求解問題。而線性代數(shù)方程組所含未知數(shù)的個數(shù)，與用來離散積分的數(shù)值方法的節(jié)點(diǎn)個數(shù)相同。在節(jié)點(diǎn)數(shù)相同的前提下，高斯數(shù)值積分方法有較高的代數(shù)精度，用它通常會得到較好的結(jié)果。對第二類Fredholm積分方程首先將積分區(qū)間[a,b]等分成n份，在每個子區(qū)間上離散方程中的積分就得到線性代數(shù)方程組。實(shí)驗(yàn)要求：分別使用如下方法，離散積分方程中的積分1.復(fù)化梯形方法；2.復(fù)化辛甫森方法；3.復(fù)化高斯方法。求解如下的積分方程，方程的準(zhǔn)確解為，并比較各算法的優(yōu)劣。程序結(jié)果：當(dāng)?shù)螖?shù)n=1時精確解復(fù)化梯形方法復(fù)化辛甫森方法復(fù)化高斯方法復(fù)化梯形方法的平均誤差err=復(fù)化辛甫森方法的平均誤差err=復(fù)化高斯方法的平均誤差err=當(dāng)?shù)螖?shù)n=5時，精確解復(fù)化梯形方法復(fù)化辛甫森方法復(fù)化高斯方法復(fù)化梯形方法的平均誤差err=復(fù)化辛甫森方法和復(fù)化高斯方法的平均誤差err=0可以看出，復(fù)化高斯方法得到的結(jié)果精度最高，復(fù)化辛普森方法比復(fù)化梯形方法的精度要高，程序：clear,clca=0;b=1;n=5;figurefun1(a,b,n);holdona=0;b=1;n=5;figurefun2(a,b,n);holdonfigurefun3(a,b,n);編制m函數(shù)：functiony=Kfun(t,s)y=2/(exp(1)-1)*exp(t);functiony=ffun(t)y=-exp(t);functiony=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;fori=1:n[x(2*i-1:2*i),w(2*i-1:2*i)]=GaussLegendre(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;fori=1:nw(2*i)=2/3*h;ifi<nw(2*i+1)=h/3;endendfunction[x,w]=fhtrapz(a,b,n)h=(b-a)/n;x=a:h:b;fori=1:n+1ifi==1||i==n+1w(i)=h/2;elsew(i)=h;endendfunction[x,w]=GaussLegendre(a,b,n)i=1:n-1;c=i./sqrt(4*i.^2-1);CM=diag(c,1)+diag(c,-1);[VL]=eig(CM);[xind]=sort(diag(L));V=V(:,ind)';w=2*V(:,1).^2;x=x*(b-a)/2+(b+a)/2;w=w*(b-a)/2;functionfun1(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);fori=1:nB(i)=ffun(x1(i));endfori=1:nforj=1:nA(i,j)=w1(j)*Kfun(x1(i),x1(j));endendA=eye(n)-A;y1=(A\B)';yN=x;fori=1:n1+1yN(i)=ffun(x(i));forj=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í)值');functionfun2(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(n,1);fori=1:nB(i)=ffun(x1(i));endfori=1:nforj=1:nA(i,j)=w1(j)*Kfun(x1(i),x1(j));endendA=eye(n)-A;y1=(A\B)';yN=x;fori=1:n1+1yN(i)=ffun(x(i));forj=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í)值');functionfun3(a,b,n)[x1,w1]=fhgauss(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);fori=1:nB(i)=ffun(x1(i));endfori=1:nforj=1:nA(i,j)=w1(j)*Kfun(x1(i),x1(j));endendA=eye(n)-A;y1=(A\B)';yN=x;fori=1:n1+1yN(i)=ffun(x(i));forj=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í)值');圖一圖二圖三實(shí)驗(yàn)三1、對常微分方程初值問題分別使用Euler顯示方法、改進(jìn)的Euler方法和經(jīng)典RK法和四階Adams預(yù)測-校正算法，求解常微分方程數(shù)值解，并與其精確解進(jìn)行作圖比較。其中多步法需要的初值由經(jīng)典RK法提供。程序：子程序functionE=Euler(fun,x0,y0,h,N)x=zeros(1,N+1);y=zeros(1,N+1);x(1)=x0;y(1)=y0;forn=1:Nx(n+1)=x(n)+h;y(n+1)=y(n)+h*feval(fun,x(n),y(n));endx1=xy1=yplot(x1,y1,'k')holdonfunctionE=Euler_modify(fun,x0,y0,h,N)x=zeros(1,N+1);y=zeros(1,N+1);x(1)=x0;y(1)=y0;forn=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);x(1)=x0;y(1)=y0;forn=1:Nx(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,,100)Euler_modify('fun',0,1,,100)[x,y]=Rk_N4('fun',0,1,,100)plot(x,y,'-*r')holdonx1=0::pi;y1=cos(x1)+sin(x1)plot(x1,y1,'')title('誤差分析');xlabel('x軸');ylabel('y軸');legend('Euler','Euler改進(jìn)法','R_K法','精確');axis([0pi]);gridon畫出圖形進(jìn)行比較：2、實(shí)驗(yàn)?zāi)康模篖orenz問題與混沌實(shí)驗(yàn)內(nèi)容：考慮著名的Lorenz方程其中s,r,b為變化區(qū)域有一定限制的實(shí)參數(shù)。該方程形式簡單，表面上看并無驚人之處，但由該方程揭示出的許多現(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的計算結(jié)果進(jìn)行對比。實(shí)驗(yàn)要求：（1）對目前取定的參數(shù)值s,r和b，選取不同的初值Y0進(jìn)行運(yùn)算，觀察計算的結(jié)果有什么特點(diǎn)解的曲線是否有界解的曲線是不是周期的或趨于某個固定點(diǎn)