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

數(shù)值分析編程作業(yè)2012年12月第二章14.考慮梯形電阻電路的設(shè)計，電路如下：電路中的各個電流{i1,i2,…,i8}須滿足以下線性方程組：這是一個三對角方程組。設(shè)V=220V，R=27，運用追趕法，求各段電路的電流量。Matlab程序如下：functionchase()%追趕法求梯形電路中各段的電流量a=input('請輸入下主對角線向量a=');b=input('請輸入主對角線向量b=');c=input('請輸入上主對角線向量c=');d=input('請輸入右端向量d=');n=input('請輸入系數(shù)矩陣維數(shù)n=');u(1)=b(1);fori=2:nl(i)=a(i)/u(i-1);u(i)=b(i)-c(i-1)*l(i);endy(1)=d(1);fori=2:ny(i)=d(i)-l(i)*y(i-1);endx(n)=y(n)/u(n);i=n-1;whilei>0x(i)=(y(i)-c(i)*x(i+1))/u(i);i=i-1;endx輸入如下:>>chase請輸入下主對角線向量a=[0,-2,-2,-2,-2,-2,-2,-2];請輸入主對角線向量b=[2,5,5,5,5,5,5,5];請輸入上主對角線向量c=[-2,-2,-2,-2,-2,-2,-2,0];請輸入方程組右端向量d=[220/27,0,0,0,0,0,0,0];請輸入系數(shù)矩陣階數(shù)n=8運行結(jié)果如下：x=8.14784.07372.03651.01750.50730.25060.11940.0477第三章14.試分別用(1)Jacobi迭代法；(2)Gauss-Seidel迭代法解線性方程組迭代初始向量?！?〕雅可比迭代法程序如下：functionjacobi()%Jacobi迭代法a=input('請輸入系數(shù)矩陣a=');b=input('請輸入右端向量b=');x0=input('請輸入初始向量x0=');n=input('請輸入系數(shù)矩陣階數(shù)n=');er=input('請輸入允許誤差er=');N=input('請輸入最大迭代次數(shù)N=');fori=1:nforj=1:nifi==jd(i,j)=a(i,j);elsed(i,j)=0;endendendm=eye(5)-d\a;%迭代矩陣g=d\b;x=m*x0+g;k=1;whilek<=N%進行迭代fori=1:5ifmax(abs(x(i)-x0(i)))>erx=m*x+g;k=k+1;elsexreturnendendcontinueendx程序執(zhí)行如下：>>jacobi請輸入系數(shù)矩陣a=[101234;19-12-3;2-173-5;32312-1;4-3-5-115]請輸入右端向量b=[12-2714-1712]'請輸入初始向量x0=[00000]'請輸入系數(shù)矩陣階數(shù)n=5請輸入允許誤差er=1.0e-6請輸入最大容許迭代次數(shù)N=60x=1.0000-2.00003.0000-2.00001.0000〔2〕高斯－賽德爾迭代法程序如下：functiongs_sdl()%gauss-seiddel迭代法a=input('請輸入系數(shù)矩陣a=');b=input('請輸入右端向量b=');x0=input('請輸入初始向量x0=');n=input('請輸入系數(shù)矩陣階數(shù)n=');er=input('請輸入允許誤差er=');N=input('請輸入最大迭代次數(shù)N=');fori=1:nforj=1:nifi<=jl(i,j)=0;elsel(i,j)=-a(i,j);endendendfori=1:nforj=1:nifi<ju(i,j)=-a(i,j);elseu(i,j)=0;endendendfori=1:nforj=1:nifi==jd(i,j)=a(i,j);elsed(i,j)=0;endendendm=(d-l)\u;%迭代矩陣g=(d-l)\b;x=m*x0+g;k=1;whilek<=Nfori=1:5ifmax(abs(x(i)-x0(i)))>erx=m*x+g;k=k+1;elsexreturnendendcontinueendx執(zhí)行結(jié)果如下：>>gs_sdl請輸入系數(shù)矩陣a=[101234;19-12-3;2-173-5;32312-1;4-3-5-115]請輸入右端向量b=[12-2714-1712]'請輸入初始向量x0=[00000]'請輸入系數(shù)矩陣階數(shù)n=5請輸入允許誤差er=1.0e-6請輸入最大容許迭代次數(shù)N=60x=1.0000-2.00003.0000-2.00001.0000第四章如下矩陣，試用冪法求按模最大的特征值與特征向量。Matlab程序代碼如下：functionmifa()A=input('請輸入系數(shù)矩陣A=');x0=input('請輸入初始列向量x0=');n=input('請輸入向量維數(shù)n=');er=input('請輸入允許誤差er=');N=input('請輸入最大容許迭代次數(shù)N=');k=1;mu=0;whilek<=Nfort=1:nifabs(x0(t))==max(abs(x0))alfa=x0(t);xb=t;%最大的x0〔i〕的下標(biāo)endendy=x0./alfa;x0=A*y;lamda=x0(xb);k=k+1;endlamda%按模最大的特征值x0%按模最大的特征值對應(yīng)的特征向量程序執(zhí)行結(jié)果如下：>>mifa請輸入系數(shù)矩陣A=[19066-8430;6630342-36;336-168147-112;30-3628291]請輸入初始列向量x0=[0001]'請輸入向量維數(shù)n=4請輸入允許誤差er=1.0e-6請輸入最大容許迭代次數(shù)N=100lamda=343.0000x0=114.3333343.0000-0.0000-171.5002第五章試編寫MATLAB函數(shù)實現(xiàn)Newton插值，要求能輸出插值多項式。對函數(shù)fx=11+4x2在區(qū)間[-5，5]上實現(xiàn)10次多項式插值。Matlab程序代碼如下：%此函數(shù)實現(xiàn)y=1/(1+4*x^2)的n次Newton插值，n由調(diào)用函數(shù)時指定%函數(shù)輸出為插值結(jié)果的系數(shù)向量〔行向量〕和插值多項式function[ty]=func5(n)x0=linspace(-5,5,n+1)';y0=1./(1.+4.*x0.^2);b=zeros(1,n+1);fori=1:n+1s=0;forj=1:it=1;fork=1:iifk~=jt=(x0(j)-x0(k))*t;end;end;s=s+y0(j)/t;end;b(i)=s;end;t=linspace(0,0,n+1);fori=1:ns=linspace(0,0,n+1);s(n+1-i:n+1)=b(i+1).*poly(x0(1:i));t=t+s;end;t(n+1)=t(n+1)+b(1);y=poly2sym(t);10次插值運行結(jié)果：[bY]=func5(10)b=Columns1through4-0.00000.00000.0027-0.0000Columns5through8-0.0514-0.00000.3920-0.0000Columns9through11-1.14330.00001.0000Y=b為插值多項式系數(shù)向量，Y為插值多項式。插值近似值：x1=linspace(-5,5,101);x=x1(2:100);y=polyval(b,x)y=Columns1through122.70033.99944.35154.09743.49262.72371.92111.17150.52740.0154-0.3571-0.5960Columns13through24-0.7159-0.7368-0.6810-0.5709-0.4278-0.2704-0.11470.02700.14580.23600.29490.3227Columns25through360.32170.29580.25040.19150.12550.0588-0.0027-0.0537-0.0900-0.1082-0.1062-0.0830Columns37through48-0.03900.02450.10520.20000.30500.41580.52800.63690.73790.82690.90020.9549Columns49through600.98861.00000.98860.95490.90020.82690.73790.63690.52800.41580.30500.2000Columns61through720.10520.0245-0.0390-0.0830-0.1062-0.1082-0.0900-0.0537-0.00270.05880.12550.1915Columns73through840.25040.29580.32170.32270.29490.23600.14580.0270-0.1147-0.2704-0.4278-0.5709Columns85through96-0.6810-0.7368-0.7159-0.5960-0.35710.01540.52741.17151.92112.72373.49264.0974Columns97through994.35153.99942.7003繪制原函數(shù)和擬合多項式的圖形代碼：plot(x,1./(1+4.*x.^2))holdallplot(x,y,'r')xlabel('X')ylabel('Y')title('Runge現(xiàn)象')gtext('原函數(shù)')gtext('十次牛頓插值多項式')繪制結(jié)果：誤差計數(shù)并繪制誤差圖：holdoffey=1./(1+4.*x.^2)-yey=Columns1through12-2.6900-3.9887-4.3403-4.0857-3.4804-2.7109-1.9077-1.1575-0.5128-0.00000.37330.6130Columns13through240.73390.75580.70100.59210.45020.29430.14010.0000-0.1169-0.2051-0.2617-0.2870Columns25through36-0.2832-0.2542-0.2053-0.1424-0.0719-0.00000.06740.12540.16960.19710.20620.1962Columns37through480.16790.12340.06600.0000-0.0691-0.1349-0.1902-0.2270-0.2379-0.2171-0.1649-0.0928Columns49through60-0.02710-0.0271-0.0928-0.1649-0.2171-0.2379-0.2270-0.1902-0.1349-0.06910.0000Columns61through720.06600.12340.16790.19620.20620.19710.16960.12540.06740.0000-0.0719-0.1424Columns73through84-0.2053-0.2542-0.2832-0.2870-0.2617-0.2051-0.11690.00000.14010.29430.45020.5921Columns85through960.70100.75580.73390.61300.37330.0000-0.5128-1.1575-1.9077-2.7109-3.4804-4.0857Columns97through99-4.3403-3.9887-2.6900plot(x,ey)xlabel('X')ylabel('ey')title('Runge現(xiàn)象誤差圖')第六章16、鋼包問題。煉鋼唱出鋼時所用的盛鋼水的鋼包，在使用過程中由于鋼液及爐渣對包襯耐火材料的侵蝕，使其容積不斷增大．經(jīng)實驗，鋼包的容積與相應(yīng)的使用次數(shù)的數(shù)據(jù)如下：使用次數(shù)x容積y使用次數(shù)x容積y106.42108.26109.58109.50109.86110.00109.93110.59110.60110.72110.90110.76111.10111.30選用雙曲線對數(shù)據(jù)進行擬合，使用最小二乘法擬合．Matlab程序如下：functiona=nihehanshu()x0=[2356791011121416171920];y0=[106.42108.26109.58109.50109.86110.00109.93110.59110.60110.72110.90110.76111.10111.30];A=zeros(2,2);B=zeros(2,1);a=zeros(2,1);x=1./x0;y=1./y0;A(1,1)=14;A(1,2)=sum(x);A(2,1)=A(1,2);A(2,2)=sum(x.^2);B(1)=sum(y);B(2)=sum(x.*y);a=A\B;y=1./(a(1)+a(2)*1./x0);subplot(1,2,2);plot(x0,y0-y,'bd-');title('擬合曲線誤差');subplot(1,2,1);plot(x0,y0,'go');holdon;x=2:0.5:20;y=1./(a(1)+a(2)*1./x);plot(x,y,'r*-');legend('散點','擬合曲線圖1/y=a(1)+a(2)*1/x');title('最小二乘法擬合曲線');求的系數(shù)為：0.00900.0008那么擬合曲線為擬合曲線圖、散點圖、誤差圖如下：第七章26.考紐螺線的形狀像鐘表的發(fā)條，也稱盤旋曲線，它在直角坐標(biāo)系中的參數(shù)方程為曲線關(guān)于原點對稱。取a=1,參數(shù)s的變化范圍[-5，5]，容許誤差限分別是10-3和10-7。選取適當(dāng)?shù)墓?jié)點個數(shù)，利用數(shù)值積分方法計算曲線上點的坐標(biāo)，并畫出曲線的圖形。程序代碼如下所示：functionhuixuan()%用梯形公式的逐次分半算法計算盤旋曲線上點的坐標(biāo)er=input('請選擇允許誤差1.0e-3或1.0e-7：');i=1;%x向量分量的下標(biāo)fors=-5:0.1:5m=1;b=s;a=0;h=(b-a)/2;fx1=cos(a^2/2);fx2=cos(b^2/2);T=h*(fx1+fx2);T0=5;whileabs(T-T0)>3*erFx=0;T0=T;fork=1:2^(m-1)%計算新增加節(jié)點處的函數(shù)值之和fx3=cos((a+(2*k-1)*h)^2/2);Fx=Fx+fx3;endT=T0/2+h*Fx;m=m+1;h=h/2;endx(i)=T;i=i+1;endj=1;%y向量分量的下標(biāo)fors=-5:0.1:5n=1;b=s;a=0;h=(b-a)/2;fy1=sin(a^2/2);fy2=sin(b^2/2);T=h*(fy1+fy2);T0=5;whileabs(T-T0)>3*erFy=0;T0=T;fork=1:2^(n-1)fy3=sin((a+(2*k-1)*h)^2/2);Fy=Fy+fy3;endT=T0/2+h*Fy;n=n+1;h=h/2;endy(j)=T;j=j+1;endplot(x,y,'k*',x,y,'k');ifer==1.0e-3title('er=1.0e-3');elsetitle('er=1.0e-7');end程序執(zhí)行結(jié)果如下：>>huixuan請選擇允許誤差1.0e-3或1.0e-7：1.0e-3>>huixuan請選擇允許誤差1.0e-3或1.0e-7：1.0e-7第八章20.求方程在附近的根，精確到。取，用簡單迭代法計算；用加快收斂的迭代格式，計算。程序及計算過程如下：建一M文件f.m存儲函數(shù)，functionf=f(x)f=exp(-x);取,用簡單迭代法計算，Matlab程序如下：function[x,i]=diedai1(x0)x=f(x0);i=1;y(i)=x;whileabs(x-x0)>10^-8i=i+1;x0=x;x=f(x);y(i)=x;end取初始值x0=0.5,輸入[x,i]=diedai1(0.5)得結(jié)果x=i=30可以看出用簡單收斂法經(jīng)過30次迭代到達精度要求。用加速收斂法的迭代格式計算，程序如下：function[x,i]=diedai2(x0)w=0.625;x=w*f(x0)+(1-w)*x0;i=1;y(i)=x;whileabs(x-x0)>10^-8i=i+1;x0=x;x=w*f(x)+(1-w)*x;y(i)=x;end同樣取x0=0.5，得x=i=5結(jié)果比擬簡單迭代法和加速迭代格式的比擬ix簡單迭代法30加速的迭代格式5可見，加速迭代格式收斂比簡單迭代格式快。第九章設(shè)有常微分方程初值問題其精確解為。選取步長使四階Adams預(yù)測-校正算法和經(jīng)典RK法均穩(wěn)定，分別用這兩種方法求解微分方程，將數(shù)值解與精確解進行比擬，輸出結(jié)果。其中多步法需要的初值由經(jīng)典RK法提供。〔1〕用經(jīng)典四階RK法求解，程序代碼如下：functionclassic_rk4()n=input('請輸入插值節(jié)點數(shù)n=');y(1)=1;f0(1)=1;%f0=cosx+sinx為精確值h=pi/n;%步長x=0:h:pi;k=2;eps=1.0e-3;fork=1:nf0(k+1)=cos(x(k))+sin(x(k));k1=-y(k)+2*cos(x(k));k2=-(y(k)+h*k1/2)+2*cos(x(k)+h/2);k3=-(y(k)+h*k2/2)+2*cos(x(k)+h/2);k4=-(y(k)+h*k3)+2*cos(x(k)+h);y(k+1)=y(k)+h/6*(k1+2*k2+2*k3+k4);endsubplot(3,1,1);plot(x,f0,'k');title('y=cosx+sinx');subplot(3,1,2);plot(x,y,'k');title('經(jīng)典四階RK法');subplot(3,1,3);T=y-f0;%計算經(jīng)典四階RK法的誤差plot(x,T,'k');title('經(jīng)典四階RK法的誤差');程序執(zhí)行結(jié)果如下：>>classic_rk4請輸入插值節(jié)點數(shù)n=3000〔2〕用四階Adams預(yù)測-校正算法求解，程序代碼如下：funct