東南大學(xué)《數(shù)值分析》上機(jī)題

東南大學(xué)《數(shù)值分析》上機(jī)題LT數(shù)值分析上機(jī)題220．（上機(jī)題）Newton迭代法（1）給定初值及容許誤差，編制Newton法解方程根的通用程序。（2）給定方程，易知其有三個根，，。1．由Newton方法的局部收斂性可知存在，當(dāng)時(shí)，Newton迭代序列收斂于根。試確定盡可能大的。2．試取若干初始值，觀察當(dāng)，，，，時(shí)Newton序列是否收斂以及收斂于哪一個根。MATLAB程序問題1clccleardx=0.5;x(1)=0.5;while(dx>1e-6)i=1;error=1;while(error>1e-8)x(i+1)=x(i)-(1/3*x(i)^3-x(i))/(x(i)^2-1);error=abs(x(i+1)-x(i));i=i+1;endif(x(i)==0)x(1)=x(1)+dx;elsedx=dx/2;x(1)=x(1)-dx;endend經(jīng)計(jì)算，最大的為0.774596問題2clcclearx2(1)=1e14;i=1;error=1;while(error>1e-8)x2(i+1)=x2(i)-(1/3*x2(i)^3-x2(i))/(x2(i)^2-1);error=abs(x2(i+1)-x2(i));i=i+1;if(i>1e4)breakendend對于不同得初始值收斂于不同的根,在（-∞，-1）內(nèi)收斂于，在（-0.774，0.774）內(nèi)收斂于0，在（1，+∞）內(nèi)收斂于，但在內(nèi)（0.774，1）和（－1，0.774）均可能收斂于和。分析：對于不同的初值，迭代序列會收斂于不同的根，所以在某個區(qū)間內(nèi)求根對于初值的選取有很大的關(guān)系。產(chǎn)生上述結(jié)果的原因是區(qū)間不滿足大范圍收斂的條件。

數(shù)值分析上機(jī)題339．（上機(jī)題）列主元三角分解法對于某電路的分析，歸結(jié)為求解線性方程組RI=V。（1）編制解n階線性方程組Ax=b的列主元三角分解法的通用程序；（2）用所編制的程序解線性方程組RI=V，并打印出解向量，保留五位有效數(shù)；（3）本編程之中，你提高了哪些編程能力？程序：clcclearA=[31,-13,0,0,0,-10,0,0,0-13,35,-9,0,-11,0,0,0,00,-9,31,-10,0,0,0,0,00,0,-10,79,-30,0,0,0,-90,0,0,-30,57,-7,0,-5,00,0,0,0,-7,47,-30,0,00,0,0,0,0,-30,41,0,00,0,0,0,-5,0,0,27,-20,0,0,-9,0,0,0,-2,29];b=[-15,27,-23,0,-20,12,-7,7,10]';[m,n]=size(A);Ap=[A,b];x=zeros(n,1);fori=1:m-1j=i;[maxa,maxi]=max(abs(Ap(i:end,j)));maxi=maxi+i-1;if(maxa~=0)mid=Ap(maxi,:);Ap(maxi,:)=Ap(i,:);Ap(i,:)=mid;fork=i:mAp(i+1:m,:)=Ap(i+1:m,:)-Ap(i+1:m,j)*(Ap(i,:)./maxa);endendendfori=linspace(m,1,m)x(i)=(Ap(i,end)-Ap(i,1:end-1)*x)/Ap(i,i);end結(jié)果：方程的解為（保留5位有效數(shù)字）：x1=-0.28923，x2=0.34544，x3=-0.71281，x4=-0.22061，x5=-0.43040，x6=0.15431，x7=-0.057823，x8=0.20105，x9=0.29023。

習(xí)題437.（上機(jī)題）3次樣條插值函數(shù)（1）編制求第一型3次樣條插值函數(shù)的通用程序;（2）已知汽車曲線型值點(diǎn)的數(shù)據(jù)如下：0123456789102.513.304.044.705.225.545.785.405.575.705.80端點(diǎn)條件為=0.8，=0.2。用所編制程序求車門的3次樣條插值函數(shù)S(x),并打印出S(i+0.5)(i=0,1,…9)。程序：（1）clcclear%%x=[0,1,2,3,4,5,6,7,8,9,10];y=[2.51,3.30,4.04,4.7,5.22,5.54,5.78,5.40,5.57,5.70,5.80];y1=0.8;yend=0.2;%%___________________________________________n=size(x,2)-1;h=x(2:end)-x(1:end-1);miu=h(1:end-1)./(h(1:end-1)+h(2:end));lamda=1-miu;f1=[y1,(y(2:end)-y(1:end-1))./h,yend];%f[xn-1,xn]f2=[f1(2:end)-f1(1:end-1)]./[h(1),h(1:end-1)+h(2:end),h(end)];%f[xn-1,xn,xn+1]A=2.*eye(n+1);A(2:end,1:end-1)=A(2:end,1:end-1)+diag([miu,1]');A(1:end-1,2:end)=A(1:end-1,2:end)+diag([1,lamda]');M=A\(6*f2');Sx=[y(1:end-1)',((y(2:end)-y(1:end-1))./h)'-((1/3*M(1:end-1)+1/6*M(2:end)).*h'),1/2*M(1:end-1),1/6*(M(2:end)-M(1:end-1))./h'];%%xx=input(’x=’);forj=2:n+1ifxx<x(j)S=Sx(j-1,:)*[1,xx-x(j-1),(xx-x(j-1))^2,(xx-x(j-1))^3]';breakendend(2)clcclear%%x=[0,1,2,3,4,5,6,7,8,9,10];y=[2.51,3.30,4.04,4.7,5.22,5.54,5.78,5.40,5.57,5.70,5.80];y1=0.8;yend=0.2;%%___________________________________________n=size(x,2)-1;h=x(2:end)-x(1:end-1);miu=h(1:end-1)./(h(1:end-1)+h(2:end));lamda=1-miu;f1=[y1,(y(2:end)-y(1:end-1))./h,yend];%f[xn-1,xn]f2=[f1(2:end)-f1(1:end-1)]./[h(1),h(1:end-1)+h(2:end),h(end)];%f[xn-1,xn,xn+1]A=2.*eye(n+1);A(2:end,1:end-1)=A(2:end,1:end-1)+diag([miu,1]');A(1:end-1,2:end)=A(1:end-1,2:end)+diag([1,lamda]');M=A\(6*f2');Sx=[y(1:end-1)',((y(2:end)-y(1:end-1))./h)'-((1/3*M(1:end-1)+1/6*M(2:end)).*h'),1/2*M(1:end-1),1/6*(M(2:end)-M(1:end-1))./h'];%%fori=0:9xx=i+0.5;forj=2:n+1ifxx<x(j)S(i+1)=Sx(j-1,:)*[1,xx-x(j-1),(xx-x(j-1))^2,(xx-x(j-1))^3]';breakendendendx∈[0,1]時(shí)；S(x)=2.51+0.8x-0.0014861x2-0.00851395x3x∈[1,2]時(shí)；S(x)=3.3+0.771486(x-1)-0.027028(x-1)2-0.00445799(x-1)3x∈[2,3]時(shí)；S(x)=4.04+0.704056(x-2)-0.0404019(x-2)2-0.0036543(x-2)3x∈[3,4]時(shí)；S(x)=4.7+0.612289(x-3)-0.0513648(x-3)2-0.0409245(x-3)3x∈[4,5]時(shí)；S(x)=5.22+0.386786(x-4)-0.174138(x-4)2+0.107352(x-4)3x∈[5,6]時(shí)；S(x)=5.54+0.360567(x-5)+0.147919(x-5)2-0.268485(x-5)3x∈[6,7]時(shí)；S(x)=5.78-0.149051(x-6)-0.657537(x-6)2+0.426588(x-6)3x∈[7,8]時(shí)；S(x)=5.4-0.184361(x-7)+0.622227(x-7)2-0.267865(x-7)3x∈[8,9]時(shí)；S(x)=5.57+0.256496(x-8)-0.181369(x-8)2+0.0548728(x-8)3x∈[9,10]時(shí)；S(x)=5.7+0.058376(x-9)-0.0167508(x-9)2+0.0583752(x-9)3S(0.5)=2.90856S(1.5)=3.67843S(2.5)=4.38147S(3.5)=4.98819S(4.5)=5.38328S(5.5)=5.7237S(6.5)=5.59441S(7.5)=5.42989S(8.5)=5.65976S(9.5)=5.7323

習(xí)題五重積分的計(jì)算23（上機(jī)題）重積分的計(jì)算題目：給定積分。取初始步長h和k，及精度。應(yīng)用復(fù)化梯形公式，采用逐次二分步長的方法，編制計(jì)算I(f)的通用程序。計(jì)算至相鄰兩次近似值之差的絕對值不超過為止。用所編程序計(jì)算積分，取。程序：clcclear%%examplef=inline('tan(x.^2+y.^2)','x','y');a=0;b=pi/3;c=0;d=pi/6;%%defineerror=1;k=1;n=1;while(error>0.5e-5)[x,y]=meshgrid(linspace(c,d,2^k+1),linspace(a,b,2^k+1));h=(b-a)/2^k;l=(d-c)/2^k;z=f(x,y);z1=z(1:end-1,1:end-1);z2=z(1:end-1,2:end);z3=z(2:end,1:end-1);z4=z(2:end,2:end);t(k)=h*l/4*(sum(sum(z1))+sum(sum(z2))+sum(sum(z3))+sum(sum(z4)));%%extrapolationif(k>=2)T(1,k-1)=4/3*t(k)-1/3*t(k-1);%T(1)error=min(error,abs(t(k)-t(k-1)));if(k>=3)T(2,k-2)=16/15*T(1,k-1)-1/15*T(1,k-2);%T(2)error=min(error,abs(T(1,k-1)-T(1,k-2)));if(k>=4)T(3,k-3)=64/63*T(1,k-2)-1/63*T(1,k-3);%T(3)error=min(error,abs(T(2,k-2)-T(2,k-3)));if(k>=5)error=min(error,abs(T(3,k-3)-T(3,k-4)));endendendendk=k+1;end計(jì)算結(jié)果：T(f)T(1)(f)T(2)(f)T(3)(f)10.51979650.3440320.3373930.33770920.38797340.3378080.3365920.3366530.35034950.3366680.3365240.33653140.34008870.3365330.33652150.33742180.33652160.3367464I(f)=0.33652二分6次

習(xí)題623．（上機(jī)題）常微分方程初值問題數(shù)值解（1）編制RK4方法的通用程序；（2）編制AB4方法的通用程序（由RK4提供初值）；（3）編制AB4-AM4預(yù)測校正方法的通用程序（由RK4提供初值）；（4）編制帶改進(jìn)的AB4-AM4預(yù)測校正方法的通用程序（由RK4提供初值）;（5）對于初值問題取步長,應(yīng)用（1）~（4）中的四種方法進(jìn)行計(jì)算，并將計(jì)算結(jié)果和精確解作比較；（6）通過本上機(jī)題，你能得到哪些結(jié)論？程序：clcclear%%Originalquestionf=inline('-x*x*y*y','x','y');y0=3;h=0.1;xstr=0;xend=1.5;x=xstr:h:xend;yx=3./(1+x.^3);n=size(x,2);%%RK4methodRK4y(1)=y0;fori=1:n-1k1=f(x(i),RK4y(i));k2=f(x(i)+h/2,RK4y(i)+h/2*k1);k3=f(x(i)+h/2,RK4y(i)+h/2*k2);k4=f(x(i)+h,RK4y(i)+h*k3);RK4y(i+1)=RK4y(i)+h/6*(k1+2*k2+2*k3+k4);end%%AB4methodAB4y(1:4)=RK4y(1:4);fori=4:n-1AB4y(i+1)=AB4y(i)+h/24*(55*f(x(i),AB4y(i))-59*f(x(i-1),AB4y(i-1))+37*f(x(i-2),AB4y(i-2))-9*f(x(i-3),AB4y(i-3)));end%%AB4-AM4predictivemethodBM4y(1:4)=RK4y(1:4);fori=4:n-1yp(i+1)=BM4y(i)+h/24*(55*f(x(i),BM4y(i))-59*f(x(i-1),BM4y(i-1))+37*f(x(i-2),BM4y(i-2))-9*f(x(i-3),BM4y(i-3)));BM4y(i+1)=BM4y(i)+h/24*(9*f(x(i+1),yp(i+1))+19*f(x(i),BM4y(i))-5*f(x(i-1),BM4y(i-1))+f(x(i-2),BM4y(i-2)));end%%ImprovedAB4-AM4predictivemethodimprBM4y(1:4)=RK4y(1:4);fori=4:n-1yP(i+1)=imprBM4y(i)+h/24*(55*f(x(i),imprBM4y(i))-59*f(x(i-1),imprBM4y(i-1))+37*f(x(i-2),imprBM4y(i-2))-9*f(x(i-3),imprBM4y(i-3)));yc(i+1)=imprBM4y(i)+h/24*(9*f(x(i+1),yP(i+1))+19*f(x(i),imprBM4y(i))-5*f(x(i-1),imprBM4y(i-1))+f(x(i-2),imprBM4y(i-2)));imprBM4y(i+1)=251/270*yc(i+1)+19/270*yP(i+1);end%%Errorerror(1:4,1:n)=abs([yx-RK4y;yx-AB4y;yx-BM4y;yx-imprBM4y]);計(jì)算結(jié)果：kx(k)y(x)RK4方法誤差A(yù)B4方法誤差A(yù)B4—AM4誤差帶改進(jìn)AB4—AM4誤差1033030303020.12.9970032.9970031.87E-072.9970031.87E-072.9970031.87E-072.9970031.87E-0730.22.976192.976193.92E-072.976193.92E-072.976193.92E-072.976193.92E-0740.32.921132.9211297.58E-072.9211297.58E-072.9211297.58E-072.9211297.58E-0750.42.8195492.8195471.61E-062.8183890.001162.8196780.000132.8195883.88E-0560.52.6666672.6666633.18E-062.6646720.0019942.6668760.0002092.6667134.62E-0570