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

數(shù)值分析上機題1,其精確值為始十七（1）編制按從大到小的順序與=￥匕+舌++，■｝，計算S,V的通用程序。（2）編制按從小到大的順序Sn=—J」；），J+白，計算Sn的通用程序。（3）按兩種順序分別計算S”，S?,S1（）6,并指出有效位數(shù)。（編制程序時用單精度）（4）通過本上機題，你明白了什么？程序代碼(matlab編程)：clccleara=single(1./([2:10A7].A2-l));SI(1)=single(0);SI(2)=1/(2A2-1);forN=3:10A2SI(N)=a(l);fori=2:N-lSI(N)=S1(N)+a(i);endendS2(1)=single(0);S2(2)=1/(2A2-1);forN=3:10A2S2(N)=a(N-l);fori=linspace(N-2,1rN-2)S2(N)=S2(N)+a(i);endendSI表示按從大到小的順序的SnS2表示按從小到大的順序的Sn計算結(jié)果從大到小的順序的值從小到大的順序的值精確值有效位數(shù)從大到小從小到大0.7400490.740050.74004965%0.7498520.74990.749944%0.7498520.7499990.74999936通過本上機題，看出按兩種不同的順序計算的結(jié)果是不相同的，按從大到小的順序計算的值與精確習題623.(上機題)常微分方程初值問題數(shù)值解(1)編制R&方法的通用程序；(2)編制AB4方法的通用程序(由RK4提供初值)；(3)編制AB4-AM4預測校正方法的通用程序(由R&提供初值)；(4)編制帶改進的AB4-AM4預測校正方法的通用程序(由R&提供初值)；(5)對于初值問題= (0<x<1.5)*0)=3取步長丸=。1,應用(1)?(4)中的四種方法進行計算，并將計算結(jié)果和精確解y(x)=3/(1+/)作比較；(6)通過本上機題，你能得到哪些結(jié)論？程序：clcclear%%Originalquestionf=inline('-x*x*y*yl'x'，'y');y0=3;h=0.1;xstr=0;xend=l.5;x=xstr:h:xend;yx=3./(l+x.A3);n=size(x,2);%%RK4methodRK4y(l)=yO;fori=l:n-1kl=f(x(i),RK4y(i));k2=f(x(i)+h/2,RK4y(i)+h/2*kl);k3=f(x(i)+h/2,RK4y(i)4-h/2*k2);k4=f(x(i)+h,RK4y(i)+h*k3);RK4y(i+l)=RK4y(i)+h/6*(kl+2*k2+2*k3+k4);end%%AB4methodAB4y(l:4)=RK4y(1:4);fori=4:n-1AB4y(i+l)=AB4y(i)+h/24*(55*f(x(i),AB4y(i))-59*f(x(i-1),AB4y(i-1))+37(x(i-2),AB4y(i-2))-9*f(x(i-3),AB4y(i-3)));end%%AB4-AM4predictivemethodBM4y(l: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)rBM4y(i-1))+37*f(x(i-2)zBM4y(i-2))-9*f(x(i-3),BM4y(i-3)));BM4y(i+1)=BM4y(i)+h/24*(9*f(x(i+l),yp(i+1))+19*f(x(i),BM4y(i))-5*f(x(i-1)zBM4y(i-1))+f(x(i-2)zBM4y(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)rimprBM4y(i-3)));yc(i+1)=imprBM4y⑴+h/24*(9*f(x(i+1),yP(i+l))+19*f(x(i),imprBM4y(i))-5*f(x(i-1),imprBM4y(i-1))+f(x(i-2),imprBM4y(i-2)));imprBM4y(i+l)=251/270*yc(i+l)+19/270*yP(i+l);end%%Errorerror(1:4f1:n)=abs([yx-RK4y;yx-AB4y;yx-BM4y;yx-imprBM4y]);計算結(jié)果：x(k)y(x) RK4方法誤差 AB4方法 誤差ab4—am4誤差帶改進ab4—am4誤差1 03 30 30303 02 0.12.997003 2.9970031.87E-07 2.9970031.87E-072.9970031.87E-072.997003 1.87E-073 0.22.97619 2.976193.92E-07 2.976193.92E-072.976193.92E-072.97619 3.92E-074 ().32.92113 2.9211297.58E-07 2.9211297.58E-072.9211297.58E-072.921129 7.58E-075 0.42.819549 2.8195471.61E-06 2.8183890.001162.8196780.000132.819588 3.88E-056 0.52.666667 2.6666633.18E-06 2.6646720.0019942.6668760.0002092.666713 4.62E-057 0.62.467105 2.46715.01E-06 2.4652030.0019032.4672520.0001472.467097 8.23E-068 0.72.233805 2.2337995.77E-06 2.2330790.0007262.2337317.35E-052.233682 ().0001229 ().81.984127 1.9841234.13E-06 1.9849510.0()08241.9837870.0()0341.983885 0.0()()24210 0.91.735107 1.7351071.16E-07 1.7370430.0019361.7346070.00051.734808 0.0002991111.5 1.5000065.81E-06 1.5021950.0021951.4995160.0004841.499732 0.00026812 1.11.287001 1.2870131.13E-05 1.2887630.0017621.2866570.0003441.286821 0.00018113 1.21.099707 1.0997221.54E-05 1.1007240.0010171.0995330.0001741.099622 8.50E-0514 1.30.93838 0.9383971.77E-05 0.93871().0003310.9383433.72E-050.938367 1.24E-0515 1.40.801282 0.80131.84E-05 0.8011350.0001470.8013274.53E-05().801311 2.93E-0516 1.50.685714 0.6857321.78E-05 0.6853350.000380.6857968.18E-050.68576 4.62E-05結(jié)論：帶改進的AB4-AM4預測校正方法比AB4-AM4預測校正方法精度更高，AB4方法精度最低,RK4方法的精度最高值有較大的誤差，而按從小到大的順序計算的值與精確值吻合。從大到小的順序計算得到的結(jié)果的有效位數(shù)少。計算機在進行數(shù)值計算時會出現(xiàn)“大數(shù)吃小數(shù)”的現(xiàn)象，導致計算結(jié)果的精度有所降低，我們在計算機中進行同號數(shù)的加法時，采用絕對值較小者先加的算法，其結(jié)果的相對誤差較小。數(shù)值分析上機題220.(上機題)Newton迭代法(1)給定初值/及容許誤差￡,編制Newton法解方程/。)=0根的通用程序。(2)給定方程/(x)=d/3—x=。，易知其有三個根耳=—JL芯=0,匕=百。.由Newton方法的局部收斂性可知存在5>0,當/￡(一夕5)時，Newton迭代序列收斂于根石。試確定盡可能大的3。.試取若干初始值，觀察當公￡(—8,7)，(T—2),(-5,6)，(氏1)，(l,oo)時Newton序列是否收斂以及收斂于哪一個根。MATLAB程序問題1clccleardx=0.5;x(1)=0.5;while(dx>le-6)i=l；error=l;while(error>le-8)x(i+1)=x(i)-(l/3*x(i)A3-x(i))/(x(i)A2-l);error=abs(x(i+1)-x(i));i=i+l;endif(x(i)==0)x(1)=x(1)+dx;elsedx=dx/2;x(1)=x(1)-dx;endend經(jīng)計算，最大的g為0.774596問題2clcclearx2(1)=lel4;i=l；error=l;while(error>le-8)x2(i+l)=x2(i)-(l/3*x2(i)"3-x2(i))/(x2(i)A2-l);error=abs(x2(i+1)-x2(i));i=i+l;if(i>le4)breakendend對于不同得初始值收斂于不同的根，/在（-8,j）內(nèi)收斂于-在（-0.774,0.774）內(nèi)收斂于0,在（1,+8）內(nèi)收斂于百，但在內(nèi)（0.774,1）和（-1,0.774）均可能收斂于—6和。分析：對于不同的初值，迭代序列會收斂于不同的根，所以在某個區(qū)間內(nèi)求根對于初值的選取有很大的關系。產(chǎn)生上述結(jié)果的原因是區(qū)間不滿足大范圍收斂的條件。數(shù)值分析上機題339.(上機題)列主元三角分解法對于某電路的分析，歸結(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,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,00f0,0r0r-5,0r0,27,-20r0r0,-9,0r0,0,-2,29];b=[-工5,27,-23,0,-20,12,-7,7,10]，[m,n]=size(A);Ap=[A,b];x=zeros(n,1);fori=l:m-l,?[maxaAmaxi]=max(abs(Ap(i:end,j)));maxi=maxi+i-l;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:mrj)*(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ù)字)：xl=-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。習題4巧0123456789102.513.304.044.705.225.545.785.405.575.705.80編制求第一型3次樣條插值函數(shù)的通用程序;已知汽車曲線型值點的數(shù)據(jù)如下：37.(上機題)編制求第一型3次樣條插值函數(shù)的通用程序;已知汽車曲線型值點的數(shù)據(jù)如下：端點條件為)，o=O.8,%。二0.2。用所編制程序求車門的3次樣條插值函數(shù)S(x),并打印出S(i+O.5)(i=0,1,-9)o程序：(1)clcclear%%x=[0,1,2,3,4,5,6,7,8,9,10];y=[2.51,3.30,4.04,4.7,5.22z5.54,5.78,5.40,5.57,5.70,5.80];yl=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=l-miu;f1=[yl,(y(2:end)-y(1:end-1))./h,yend];%f[xn-1,xn]f2=[fl(2:end)-fl(1:end-1)]./[h(l),h(l: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]1);A(1:end-1z2:end)=A(1:end-1,2:end)+diag([1,lamda]');M=A\(6*f21);Sx=[y(1:end-1)1,((y(2:end)-y(1:end-1))./h),-((1/3*M(1:end-1)+l/6*M(2:end)).*hf),1/2*M(1:end-l),1/6*(M(2:end)-M(l:end-l))./h1]；%%xx=input('x=');forj=2:n+lifxx<x(j)S=Sx(j-1,:)*[1,xx-x(j-1)z(xx-x(j-1))A2,(xx-x(j-1))^3]1;breakendendclcclcclear%%x=[0,1,2,3,4,5,6,7,8,9,10];y=[2.51,3.30,4.04,4.7,5.22z5.54,5.78,5.40,5.57,5.70,5.80];yl=0.8;yend=0.2;%%n=size(xz2)-1;h=x(2:end)-x(1:end-1);miu=h(1:end-1)./(h(1:end-1)+h(2:end));lamda=l-miu;f1=[yl,(y(2:end)-y(1:end-1))./h,yend];%f[xn-1,xn]f2=[fl(2:end)-fl(1:end-1)]./[h(l),h(l: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-1z2:end)=A(1:end-1z2:end)+diag([1,lamda]');M=A\(6*f21);Sx=[y(1:end-1)1z((y(2:end)-y(1:end-1))./h)1-((1/3*M(1:end-1)+l/6*M(2:end)).)，1/2*M(1:end-l),1/6*(M(2:end)-M(l:end-1))./h1]；%%fori=0:9xx=i+0.5;forj=2:n+1ifxx<x(j)S(i4-l)=Sx(j-l,:)*[l,xx-x(j-l),(xx-x(j-l))A2Z(xx-x(j-l))A3]f;breakendendendxe[0,1]時;2一0.00851395x3x￡[1,2]時；S(x)=3.3+0.771486(x-1)-0.027028(x-1)2-0.00445799(x-1)3xe[2,3]時；S(x)=4.04+0.704056(x-2)-0.0404019(x-2)2-0.0036543(x-2)3xe[3,4]時；S(x)=4.7+0.612289(x-3)-0.0513648(x-3)2-0.0409245(x-3)3x￡[4,5]時;S(x)=5.22+0.386786(x-4)-0.174138(x-4)2+0.107352(x-4)3x￡[5,6]時；S(x)=5.54+0.360567(x-5)+0.147919(x-5)2-0.268485(x-5)3xe[6,7]時；S(x)=5.78-0.14905l(x-6)-0.657537(x-6)2+0.426588(x-6)3xe[7,8]時；S(x)=5.4-0.18436l(x-7)+0.622227(x-7)2-0.267865(x-7)3xe[8,9]時;S(x)=5.57+0.256496(x-8)-0.181369(x-8)2+0.0548728(x-8)3xe[9,l0]時；S(x)=5.7+0.058376(x-9)-0.0167508(x-9)2+0.0583752(x-9)3S(0.5>2.90856 S(l.5)=3.67843 S(2.5)=4.38147S(3.5)=4.98819 S(4.5)=5.38328 S(5.5)=5.7237S(6.5)=5.59441 S(7.5)=5.42989 S(8.5)=5.65976S(9.5)=5.7323習題五重積分的計算23(上機題)重積分的計算題目：給定積分/(/)=1(，/(x,y)dx)dy。取初始步長h和k,及精度￡。應用復化梯形公式，采用逐次二分步長的方法，編制計算I⑴的通用程序。計算至相鄰兩次近似值之差的絕對值不超過￡為止。1)用所編程序計算積分/(7)=『(『火，+p2)公)力，取￡=0.5*10-5。程序：clcclear%%examplef=inline(!tan(x.人2+y.人2)，,，x，,，y。;a=0;b=pi/3;c=0;d=pi/6;%%define