數(shù)值分析報(bào)告

數(shù)值分析實(shí)驗(yàn)報(bào)告學(xué)院：土木與建筑學(xué)院姓名：馬X學(xué)號(hào)：10XXXXXXX專業(yè)：xxxxxxxxxxxx班第二章插值法一實(shí)驗(yàn)?zāi)康?.了解拉格朗日插值法，三次樣條插值，牛頓插值法的基本原理和方法。2.通過(guò)實(shí)例掌握用MATLAB求插值的方法。3.編程實(shí)現(xiàn)拉格朗日插值，三次樣條插值，牛頓插值。二實(shí)驗(yàn)內(nèi)容1已知函數(shù)在下列各點(diǎn)的值為xi0.20.40.60.81.0f(xi)0.980.920.810.640.38試用4次牛頓插值多項(xiàng)式P4（x）及三次樣條函數(shù)S（x）（自然邊界條件）對(duì)數(shù)據(jù)進(jìn)行插值。用圖給出{（xi,yi）,xi=0.2+0.8i,i=0,1,11,10},P4(x)及S（x）.2．在區(qū)間[-1,1]上分別取n=10，20用兩組等距節(jié)點(diǎn)對(duì)龍格函數(shù)f(x)=1/(1+25x2)作多項(xiàng)式插值及三次樣條插值，對(duì)每個(gè)n值，分別畫出插值函數(shù)及f(x)的圖形。三．實(shí)驗(yàn)步驟1開啟軟件平臺(tái)---MATLAB，打開MATLAB編輯窗口。2.根據(jù)算法編制程序并運(yùn)行(1).編制牛頓插值程序：程序代碼：disp('牛頓插值');x=0.2:0.2:1.0;y=[0.980.920.810.640.38];plot(x,y,'r');holdon;xi=[0.20.281.081.0];n=length(x);m=length(y);ifm~=nerror('向量x和向量y的長(zhǎng)度必須一致')endw=ones(1,length(xi));a(1)=y(1);yy=y(1);fork=2:n;a(k)=y(k);w=w.*(xx-x(k-1));fori=1:k-1;a(k)=(a(i)-a(k))/(x(i)-x(k));endyi=yi+w*a(k);endyiplot(xi,yi,'k');（以上是牛頓插值多項(xiàng)式）運(yùn)行>>x=0.2:0.2:1.0;>>y=[0.980.920.810.640.38];holdon;xx=[0.20.281.081.0];yy=spline(x,y,xx);plot(x,y,'o',xx,yy);pp=csape(x,y,'variational',[00]);yy=ppval(pp,xx);fnplt(pp);yy(以上是三次樣條插值)（2）編制程序代碼：拉格朗日插值functionyy=nalagr(x,y,xx)m=length(x);n=length(y);ifm~=n,error('x和y的向量長(zhǎng)度必須相等');ends=0;fori=1:nt=ones(1,length(xx));forj=1:nifj~=i,t=t.*(xx-x(j))/(x(i)-x(j));endends=s+t*y(i);endyy=s;運(yùn)行>>x=-1:1y=1./(1+x.^2);xx=-1:0.1:1;yy=nalagr(x,y,xx);plot(x,y,'o',xx,yy);>>x=-1:1y=1./(1+x.^2);xx=-1:0.2:1;yy=nalagr(x,y,xx);plot(x,y,'o',xx,yy);三次樣條插值：n=10運(yùn)行>>x=-1:1;y=(1+25*x.^2).^(-1);xx=-1:0.2:1;yy=spline(x,y,xx);plot(x,y,'o',xx,yy);n=20運(yùn)行>>x=-1:1;y=(1+25*x.^2).^(-1);xx=-1:0.1:1;yy=spline(x,y,xx);plot(x,y,'o',xx,yy);第三章函數(shù)逼近與快速傅里葉變換一實(shí)驗(yàn)?zāi)康?.了解最小二乘擬合的基本原理和方法；2.掌握用MATLAB作最小二乘多項(xiàng)式擬合和曲線擬合的方法；3.通過(guò)實(shí)例學(xué)習(xí)如何用擬合方法解決實(shí)際問(wèn)題，注意與插值方法的區(qū)別。4.了解各種參數(shù)辨識(shí)的原理和方法；5.通過(guò)范例展現(xiàn)由機(jī)理分析確定模型結(jié)構(gòu)，擬合方法辨識(shí)參數(shù)，誤差分析等求解實(shí)際問(wèn)題的過(guò)程。二實(shí)驗(yàn)內(nèi)容1.對(duì)于給函數(shù)f(x)=1/(1+25x2)在區(qū)間[-1,1]上取xi=-1+0.2i(i=0,1,………10),試求3次曲線擬合，試畫出擬合曲線并打印出方程，與第2章計(jì)算實(shí)習(xí)題2的結(jié)果比較。2.由實(shí)驗(yàn)給出數(shù)據(jù)表x0.00.10.20.30.50.81.0y1.00.410.500.610.912.022.46試求3次，4次多項(xiàng)式的曲線擬合，再根據(jù)曲線形狀，求另外函數(shù)的擬合曲線，用圖示數(shù)據(jù)曲線及相應(yīng)的三種擬合曲線。三實(shí)驗(yàn)步驟1開啟軟件平臺(tái)----MATLAB，打開MATLAB編輯窗口。2根據(jù)算法編制程序，并運(yùn)行，觀察結(jié)果（1）編制程序曲線擬合程序；functionp=nafit(x,y,m)A=zeros(m+1,m+1);fori=0:mforj=0:mA(i+1,j+1)=sum(x.^(i+j));endb(i+1)=sum(x.^i.*y);enda=A\b';p=fliplr(a');plot(x,y);holdon運(yùn)行>>x=-1:0.2:1;y=(1+25*x.^2).^(-1);nafit(x,y,3);ans=0.0000-0.5752-0.00000.4841（2）直接引用（1）的程序。運(yùn)行：>>x=[0.0,0.1,0.2,0.3,0.5,0.8,1.0];y=[1.0,0.41,0.50,0.61,0.91,2.02,2.46];nafit(x,y,3)ans=-6.622112.8147-4.65910.9266>>x=[0.0,0.1,0.2,0.3,0.5,0.8,1.0];y=[1.0,0.41,0.50,0.61,0.91,2.02,2.46];nafit(x,y,4)ans=2.8853-12.334816.2747-5.29870.9427>>x=[0.0,0.1,0.2,0.3,0.5,0.8,1.0];y=[1.0,0.41,0.50,0.61,0.91,2.02,2.46];p=polyfit(x,y,5);xi=-0.0:0.05:1.0;p;yi=polyval(p,xi);plot(x,y,'o',xi,yi,'k');pp=-79.3261195.4554-172.710469.0498-11.00440.9955第四章數(shù)值積分與數(shù)值微分一實(shí)驗(yàn)?zāi)康?.了解數(shù)值積分的基本原理和方法；2.了解書上介紹的幾種數(shù)值積分的不同原理和方法；3.掌握用MATLAB求積分方法：復(fù)合梯形法，復(fù)合辛普森法，龍貝格法及自適應(yīng)辛普森法。4.通過(guò)實(shí)例學(xué)習(xí)如何用以上幾種方法求積分。二實(shí)驗(yàn)內(nèi)容：用不同數(shù)值方法計(jì)算積分。（1）取不同的步長(zhǎng)h，分別用復(fù)合梯形及復(fù)合辛普森求積計(jì)算積分，給出誤差中關(guān)于h的函數(shù)，并與積分精確值比較兩個(gè)公式的精度，是否存在一個(gè)最小的h，使得精度不能再被改善？（2）用龍貝格求積計(jì)算完成問(wèn)題（1）。（3）用自適應(yīng)辛普森積分，試其精度達(dá)到10-4.三實(shí)驗(yàn)步驟1開啟軟件平臺(tái)MATLAB，打開MATLAB編輯窗口。2根據(jù)不同積分算法編制程序，并運(yùn)行，觀察結(jié)果。（1）復(fù)合梯形求積functiont=natrapz(fname,a,b,n)h=(b-a)/n;fa=feval(fname,a);fb=feval(fname,b);f=feval(fname,a+h:h:b-h+0.001*h);t=h*(0.5*(fa+fb)+sum(f));運(yùn)行>>formatlong;natrapz(inline('sqrt(x).*log(x)'),eps,1,16),formatshort;ans=-0.42947460162951復(fù)合辛普森functiong=nagsint(fname,a,b,n,m)switchmcase1t=0,A=2;case2t=[-1/sqrt(3),1/sqrt(3)];A=[1,1];case3t=[-sqrt(0.6),0,sqrt(0.6)];A=[5/9,8/9,5/9];case4t=[-0.861136-0.3399810.3399810.861136];A=[0.3478550.6521450.6521450.347855];case5t=[-0.906180-0.53846900.5384690.906180];A=[0.2369270.4786290.5688890.4786290.236927];otherwiseerror;endx=linspace(a,b,n+1);g=0;fori=1:ng=g+gsint(fname,x(i),x(i+1),A,t);endfunctiong=gsint(fname,a,b,A,t)g=(b-a)/2*sum(A.*feval(fname,(b-a)/2*t+(a+b)/2));運(yùn)行>>formatlong;nagsint(inline('sqrt(x).*log(x)'),eps,1,2,3),formatshort;ans=-0.44797782008964（2）龍貝格求積程序functiont=naromberg(fname,a,b,e)ifnargin<2,e=1e-4;end;i=1;j=1;h=b-a;T(i,1)=h/2*(feval(fname,a)+feval(fname,b));T(i+1,1)=T(i,1)/2+sum(feval(fname,a+h/2:h:b-h/2+0.001*h))*h/2;T(i+1,j+1)=4^j*T(i+1,j)/(4^j-1)-T(i,j)/(4^j-1);whileabs(T(i+1,i+1)-T(i,i))>ei=i+1;h=h/2;T(i+1,1)=T(i,1)/2+sum(feval(fname,a+h/2:h:b-h/2+...0.001*h))*h/2;forj=1:iT(i+1,j+1)=4^j*T(i+1,j)/(4^j-1)-T(i,j)/(4^j-1);endendT；t=T(i+1,j+1);運(yùn)行：>>formatlong;naromberg(inline('sqrt(x).*log(x)'),eps,1,0.5e-6),formatshort;ans=-0.44444428266070（3）自適應(yīng)辛普森積分程序：function[q,srm,err]=naadapt(fname,a,b,e)srm=zeros(30,6);it=0;done=1;m=1;sr=simpson(fname,a,b,e);srm(1,1:6)=sr;state=it;while(state==it)n=m;forj=n:-1:1p=j;sr0=srm(p,:);err=sr0(5);e=sr0(6);ife<=errstate=done;a=sr0(1);b=sr0(2);err=sr0(5);e=sr0(6);c=0.5*(a+b);e2=0.5*e;sr1=simpson(fname,a,c,e2);sr2=simpson(fname,c,b,e2);err=abs(sr0(3)-sr1(3)-sr2(3))/10;iferr<esrm(p,:)=sr0;srm(p,4)=sr1(3)+sr2(3);srm(p,5)=err;elsesrm(p+1:m+1,:)=srm(p:m,:);m=m+1;srm(p,:)=sr1;srm(p+1,:)=sr2;state=it;endendendendq=sum(srm(:,4));err=sum(abs(srm(:,5)));srm=srm(1:m,:);functionz=simpson(fname,a,b,e)h=b-a;c=0.5*(a+b);f=feval(fname,[acb]);s=h*(f(1)+4*f(2)+f(3))/6;s2=s;err=e;z=[abss2erre];運(yùn)行：>>[q,srm,err]=naadapt(inline('sqrt(x).*log(x)'),eps,1,1e-4)q=-0.4444srm=0.00000.0000-0.0000-0.00000.00000.00000.00000.0000-0.0000-0.00000.00000.00000.00000.0000-0.0000-0.00000.00000.00000.00000.0001-0.0000-0.00000.00000.00000.00010.0001-0.0000-0.00000.00000.00000.00010.0002-0.0000-0.00000.00000.00000.00020.0005-0.0000-0.00000.00000.00000.00050.0010-0.0001-0.00010.00000.00000.00100.0020-0.0002-0.00020.00000.00000.00200.0039-0.0006-0.00060.00000.00000.00390.0078-0.0015-0.00150.00000.00000.00780.0156-0.0037-0.00370.00000.00000.01560.0313-0.0089-0.00890.00000.00000.03130.0625-0.0206-0.02060.00000.00000.06250.1250-0.0451-0.04510.00000.00000.12500.2500-0.0902-0.09020.00000.00000.25000.5000-0.1494-0.14940.00000.00000.50001.0000-0.1239-0.12390.00000.0001err=7.7219e-006第五章解線性方程組的直接解法一實(shí)驗(yàn)?zāi)康?.了解矩陣范數(shù)的定義和常見(jiàn)的幾種矩陣范數(shù)；2.掌握用直接法解線性方程組的常見(jiàn)幾種方法；3.通過(guò)實(shí)例學(xué)習(xí)如何使用直接法中的LU分解及列主高斯消去法求解線性方程組。二實(shí)驗(yàn)內(nèi)容用LU分解及列主高斯消去法解線性方程組=輸出Ax=b中系數(shù)A=LU分解的矩陣L及U，解向量x及detA;列主元法的行交換次序，解向量x及detA；比較兩種方法所得的結(jié)果。2用列主元高斯消去法解線性方程組Ax=b，（1）（2）分別輸出A，b，detA,解向量x，（1）中A的條件數(shù)。分析比較（1），（2）的計(jì)算結(jié)果。線性方程組Ax=b的A及b為，，則解x=（1，1，1，1）T。用MATLAB內(nèi)部函數(shù)求detA及A的所有特征值和cond(A)2.若令，求解（A+﹠A）（x+﹠x）=b，輸出向量﹠x和。從理論結(jié)果和實(shí)際計(jì)算兩方面分析線性方程組Ax=b解得相對(duì)誤差2/及A的相對(duì)誤差和/的關(guān)系。三實(shí)驗(yàn)步驟1.開啟軟件平臺(tái)MATLAB，打開MATLAB編輯窗口。2.根據(jù)計(jì)算方法及原理編輯程序，運(yùn)行，并觀查結(jié)果。（1）第一題LU分解程序：function[l,u]=nalu(a)n=length(a);u=zeros(n,n);l=eye(n,n);u(1,:)=a(1,:);l(2:n,1)=a(2:n,1)/u(1,1);fork=2:nu(k,k:n)=a(k,k:n)-l(k,1:k-1)*u(1:k-1,k:n);l(k+1:n,k)=(a(k+1:n,k)-l(k+1:n,1:k-1)*u(1:k-1,k))/u(k,k);end運(yùn)行：>>a=[10-701;-32.09999962;5-15-1;2102];[l,u]=nalu(a)l=1.0e+006*0.0000000-0.00000.0000000.0000-2.50000.000000.0000-2.40000.00000.0000u=1.0e+007*0.0000-0.000000.00000-0.00000.00000.0000001.50000.57500000.0000>>x=(u^(-1))*(l^(-1))*[8;5.900001;5;1]x=0.0000-1.00001.00001.0000>>d=det(a)d=-762.0001列主高斯消元法程序functionx=nagauss2(a,b,flag)ifnargin<3,flag=0;endn=length(b);a=[a,b];fork=1:(n-1)[ap,p]=max(abs(a(k:n,k)));p=p+k-1;ifp>k,t=a(k,:);a(k,:)=a(p,:);a(p,:)=t;enda((k+1):n,(k+1):(n+1))=a((k+1):n,(k+1):(n+1))...-a((k+1):n,k)/a(k,k)*a(k,(k+1):(n+1));a((k+1):n,k)=zeros(n-k,1);ifflag==0,a;endendx=zeros(n,1);x(n)=a(n,n+1)/a(n,n);fork=n-1:-1:1x(k,:)=(a(k,n+1)-a(k,(k+1):n)*x((k+1):n))/a(k,k);end運(yùn)行：>>a=[10-701;-32.09999962;5-15-1;2102];b=[8;5.900001;5;1];x=nagauss2(a,b)x=0.0000-1.00001.00001.0000>>d=det(a)d=-762.0001（2）第2題第一小題引用題1題的列主高斯消元程序并運(yùn)行。>>a=[3.016.031.99;1.274.16-1.23;0.987-4.819.34];b=[111]';x=nagauss2(a,b)x=1.0e+003*1.5926-0.6319-0.4936>>det(a)ans=-0.0305第二小題引用第1題的列主高斯消元程序并運(yùn)行。>>a=[3.006.031.99;1.274.16-1.23;0.990-4.819.34];b=[111]';x=nagauss2(a,b)x=119.5273-47.1426-36.8403>>det(a)ans=-0.4070第6章解線性方程組的迭代解一實(shí)驗(yàn)?zāi)康?.了解求解線性方程的解的常見(jiàn)方法；2.掌握用雅可比迭代法和SOR迭代法解線性方程組；3.通過(guò)實(shí)例學(xué)習(xí)用MATLAB使用雅可比迭代法和SOR迭代法求解線性方程。4.掌握判斷雅可比迭代法和SOR迭代法收斂的條件，和收斂階的計(jì)算二實(shí)驗(yàn)內(nèi)容給出線性方程組，其中系數(shù)矩陣為希爾伯特矩陣：.假設(shè)若取，分別用雅可比迭代法及SOR迭代（）求解.比較計(jì)算結(jié)果。三實(shí)驗(yàn)步驟開啟軟件平臺(tái)MATLAB，打開matlab編輯窗口。編寫程序代碼，運(yùn)行并觀查代碼若，在命令窗口輸入：>>H=hilb(6)%：H=hilb(6)的作用是生成6階希爾伯特矩陣。運(yùn)行結(jié)果為：H=1.00000.50000.33330.25000.20000.16670.50000.33330.25000.20000.16670.14290.33330.25000.20000.16670.14290.12500.25000.20000.16670.14290.12500.11110.20000.16670.14290.12500.11110.10000.16670.14290.12500.11110.10000.0909為了求，在命令窗口中輸入：>>x=[1;1;1;1;1;1];>>b=H*x運(yùn)行結(jié)果為：b=2.45001.59291.21790.99560.84560.7365采用雅可比迭代法，程序?yàn)椋篺unction[y,n]=jacobi(A,b,x0,eps)ifnargin==3eps=1.0e-6;elseifnargin<3errorreturnendD=diag(diag(A));%求A的對(duì)角矩陣L=-tril(A,-1);%求A的下三角陣U=-triu(A,1);%求A的上三角陣B=D\(L+U);f=D\b;y=B*x0+f;n=1;%迭代次數(shù)whilenorm(y-x0)>=epsx0=y;y=B*x0+f;n=n+1;end求解方程組時(shí)，在命令窗口中輸入：y=InfInfNaNNaNNaNNaNn=487可見(jiàn)雅可比迭代法不收斂。采用SOR迭代，程序?yàn)椋篺unctionx=nasor(A,b,omega,x0,e,N)%：x=nasor(A,b,omega,x0,e,N)A為系數(shù)矩陣，b為右端向量，x返回解向量，%x0為初值向量（默認(rèn)原點(diǎn)），e為精度（默認(rèn)1e-4），設(shè)置迭代次數(shù)上限以防發(fā)散%（默認(rèn)500），omega是松弛因子，默認(rèn)為1.5n=length(b);ifnargin<6,N=500;endifnargin<5,e=1e-4;endifnargin<4,x0=zeros(n,1);endifnargin<3,omega=1.5;endx=x0;x0=x+2*e;k=0;L=tril(A,-1);U=triu(A,1);whilenorm(x0-x,inf)>e&k<N,k=k+1,x0=x;fori=1:nx1(i)=(b(i)-L(i,1:i-1)*x(1:i-1,1)-U(i,i+1:n)*x0(i+1:n,1))/A(i,i);x(i)=(1-omega)*x0(i)+omega*x1(i);enddisp(x')endifk==N,warning('已達(dá)迭代次數(shù)上限');end求解方程組時(shí)，當(dāng)時(shí)在命令窗口中輸入：>>nasor(H,b,1,[0,0,0,0,0,0]’,1e-6),formatshort當(dāng)時(shí)，運(yùn)行結(jié)果為：......k=500Columns1through51.001692280872260.995319364990660.973515481781651.049050118497851.02646662676324Column60.95160554693359Warning:已達(dá)到迭代次數(shù)上限>Innasorat17ans=1.001692280872260.995319364990660.973515481781651.049050118497851.026466626763240.95160554693359當(dāng)時(shí)，運(yùn)用雅可比迭代法均不收斂。當(dāng)，，迭代初值為[0,0,0,0,0,0,0,0]T，用SOR方法求得方程組的解為：當(dāng)，，迭代初值為[0,0,0,0,0,0,0,0,0,0]T用SOR方法求得方程組的解為：第七章非線性方程與方程組的數(shù)值解法一實(shí)驗(yàn)?zāi)康?.了解求解非線性方程的解的常見(jiàn)方法；2.掌握用不動(dòng)點(diǎn)迭代法和牛頓迭代法解線性方程組；3.通過(guò)實(shí)例學(xué)習(xí)用MATLAB使用不動(dòng)點(diǎn)迭代法和牛頓迭代法及斯特芬森加速迭代法分解求解非線性方程。4.掌握判斷不動(dòng)點(diǎn)迭代法和牛頓迭代法收斂的條件，和收斂階的計(jì)算；二實(shí)驗(yàn)內(nèi)容求下列方程的實(shí)根：x2-3x+2-ex=0;(2)x3+2x2+10x-20=0.三實(shí)驗(yàn)要求（1）設(shè)計(jì)一種不動(dòng)點(diǎn)迭代法，要使迭代序列收斂，然后再用斯特芬森加速迭代，計(jì)算到｜xk-xk-1｜<10-8為止。（2）用牛頓迭代，同樣計(jì)算到｜xk-xk-1｜<10-8輸出迭代初值及各次迭代值和迭代次數(shù)k，比較方法的優(yōu)劣。四實(shí)驗(yàn)步驟開啟軟件平臺(tái)MATLAB，打開MATLAB編輯窗口。根據(jù)迭代法的計(jì)算方法和原理，編寫不動(dòng)點(diǎn)迭代程序和牛頓迭代程序及斯特芬森加速迭代。（1）第1小題>>symsxf=x^2-3*x+2-exp(x);df=diff(f,x);eps=1e-8;x0=6;fprintfx0=6;cnt=0;MAXCNT=200;whilecnt<MAXCNTx1=x0-subs(f,x,x0)/subs(df,x,x0)if(abs(x1-x0)<eps)break;endx0=x1;cnt=cnt+1;endifcnt==MAXCNTdispelsevpa(x1,8)end運(yùn)行x0=6x1=5.0279x1=4.0632x1=3.0857x1=2.0372x1=0.8802x1=0.2565x1=0.2575x1=0.2575x1=0.2575ans=.25753029（2）第2題>>symsxf=x^3+2*x^2+10*x-20;df=diff(f,x);eps=1e-8;x0=6;fprintfx0=6;cnt=0;MAXCNT=200;whilecnt<MAXCNTx1=x0-subs(f,x,x0)/subs(df,x,x0)if(abs(x1-x0)<eps)break;endx0=x1;cnt=cnt+1;endifcnt==MAXCNTdispelsevpa(x1,8)end運(yùn)行x0=6x1=3.6901x1=2.2516x1=1.5481x1=1.3777x1=1.3688x1=1.3688x1=1.3688第九章常微分方程初值問(wèn)題數(shù)值解法一實(shí)驗(yàn)?zāi)康氖煜じ鞣N初值問(wèn)題的算法，用MATLAB編出算法程序；明確各種算法的精度與所選步長(zhǎng)有密切關(guān)系；通過(guò)計(jì)算更加了解各種算法的優(yōu)越性。二實(shí)驗(yàn)內(nèi)容y’=1/x2-y/x,1≤x≤2,y(1)=1;y’=-50y+50x2+2x,0≤x≤1.考慮常微分方程組初值問(wèn)題其中，三實(shí)驗(yàn)要求用改進(jìn)歐拉法（h=0.05）及經(jīng)典四階R-K法（h=0.1）求（1）的數(shù)值解，并打印x=1+0.1i（i=0,1,…….10）的值。用經(jīng)典四階R-K方法解（2），步長(zhǎng)分別取h=0.1,0.025,0.01,計(jì)算并打印x=0.1i(i=0,1,…10)各點(diǎn)的值，與準(zhǔn)確解y(x)=1/3e-50x+x2比較。要求用四階R-K方法及梯形法計(jì)算（可以直接用數(shù)學(xué)庫(kù)的軟件），根據(jù)計(jì)算結(jié)果畫出函數(shù)的圖形。四實(shí)驗(yàn)步驟開啟軟件平臺(tái)MATLAB，打開編輯窗口。根據(jù)求數(shù)值解得計(jì)算方法和原理，編寫改進(jìn)歐拉法和經(jīng)典四階R-K法的程序代碼。運(yùn)行并觀查結(jié)果。（a）改進(jìn)歐拉法程序function[x,y]=naeuler2(dyfun,xspan,y0,h)x=xspan(1):h:xspan(2);y(1)=y0;forn=1:length(x)-1k1=feval(dyfun,x(n),y(n));y(n+1)=y(n)+h*k1;k2=feval(dyfun,x(n+1),y(n+1));y(n+1)=y(n)+h*(k1+k2)/2;endx=x';y=y';運(yùn)行：>>clear;dyfun=inline('1/(x^2)-y/x');[x,y]=naeuler2(dyfun,[1,2],1,0.05);[x,y]ans=1.00001.00001.05000.99891.10000.99581.15000.99111.20000.98531.25000.97861.30000.97111.35000.96311.40000.95471.45000.94601.50000.93711.55000.92801.60000.91881.65000.90961.70000.90041.75000.89131.80000.88221.85000.87321.90000.86421.95000.85542.00000.8467經(jīng)典四階程序代碼：function[x,y]=nark4(dyfun,xspan,y0,h)x=xspan(1):h:xspan(2);y(1)=y0;forn=1:length(x)-1k1=feval(dyfun,x(n),y(n));k2=feval(dyfun,x(n)+h/2,y(n)+h/2*k1);k3=feval(dyfun,x(n)+h/2,y(n)+h/2*k2);k4=feval(dyfun,x(n+1),y(n)+h*k3);y(n+1)=y(n)+h*(k1+2*k2+2*k3+k4)/6;endx=x';y=y';運(yùn)行：>>clear;dyfun=inline('1/(x^2)-y/x');[x,y]=nark4(dyfun,[1,2],1,0.1);[x,y]ans=1.00001.00001.10000.99571.20000.98531.30000.97101.40000.95461.50000.93701.60000.91881.70000.90041.80000.88211.90000.86412.00000.8466（b）h=0.1時(shí)，引用（1）編寫的R-K程序代碼。>>clear;dyfun=inline('-50*y+50*x^2+2*x');[x,y]=nark4(dyfun,[0,1],1/3,0.1);[x,y]ans=1.0e+010*00.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00000.00020.00000.00310.00000.04180.00000.57330.00007.8594h=0.025時(shí)>>clear;dyfun=inline('-50*y+50*x^2+2*x');[x,y]=nark4(dyfun,[0,1],1/3,0.025);[x,y]ans=00.33330.02500.10310.05000.03400.07500.01530.10000.01300.12500.01660.15000.02280.17500.03070.20000.04010.22500.05070.25000.06250.27500.07570.30000.09000.32500.10570.35000.12250.37500.14070.40000.16000.42500.18070.45000.20250.47500.22570.50000.25000.52500.27570.55000.30250.57500.33070.60000.36000.62500.39070.65000.42250.67500.45570.70000.49000.72500.52570.75000.56250.77500.60070.80000.64000.82500.68070.85000.72250.87500.76570.90000.81000.92500.85570.95000.90250.97500.95071.00001.0000h=0.01時(shí)，引用（1）編寫的R-K程序代碼。clear;dyfun=inline('-50*y+50*x^2+2*x');[x,y]=nark4(dyfun,[0,1],1/3,0.01);[x,y]ans=00.33330.01000.20240.02000.12310.03000.07540.04000.04680.05000.02990.06000.02020.07000.01500.08000.01250.09000.01180.10000.01230.11000.01350.12000.01520.13000.01740.14000.01990.15000.02270.16000.02570.17000.02900.18000.03240.19000.03610.20000.04000.21000.04410.22000.04840.23000.05290.24000.05760.25000.06250.26000.06760.27000.07290.28000.07840.29000.08410.30000.09000.31000.09610.32000.10240.33000.10890.34000.11560.35000.1225