數(shù)值方法計(jì)算實(shí)習(xí)題1

信計(jì)091龔立麗200900901004數(shù)值方法計(jì)算實(shí)習(xí)題要求：1、用Matlab語言或你熟悉的其他算法語言編寫程序，使之盡可能具有通用性；2、根據(jù)上機(jī)計(jì)算實(shí)踐，對(duì)所使用的數(shù)值方法的特點(diǎn)、性質(zhì)、有效性、誤差和收斂性等方面進(jìn)行必要的討論和分析；3、完成計(jì)算后寫出實(shí)驗(yàn)報(bào)告，內(nèi)容包括：課題名稱、解決的問題、采用的數(shù)值方法、算法程序、數(shù)值結(jié)果、對(duì)實(shí)驗(yàn)結(jié)果的討論和分析等；4、特別說明：嚴(yán)禁抄襲，否那么一經(jīng)發(fā)現(xiàn)，所有雷同實(shí)驗(yàn)報(bào)告最多評(píng)為及格。一、下表給出了飛行中鴨子的上部形狀的節(jié)點(diǎn)數(shù)據(jù)，試用三次樣條插值函數(shù)〔自然邊界條件〕和20次Lagrange插值多項(xiàng)式對(duì)數(shù)據(jù)進(jìn)行插值。用圖示出給定的數(shù)據(jù)，以及和。0.91.31.92.12.63.03.94.44.75.06.01.31.51.852.12.62.72.42.152.052.12.257.08.09.210.511.311.61212.613.013.32.32.251.951.40.90.70.60.50.40.25解：>>x=[0.91.31.92.12.63.03.94.44.75.06.07.08.09.210.511.311.61212.613.013.3];>>y=[1.31.51.852.12.62.72.42.152.052.12.252.32.251.951.40.90.70.60.50.40.25];〔1〕三次樣條插值法在MATLAB中編寫m文件function[f,f0]=scyt(x,y,y2_1,y2_N,x0)%y2_1和y2_N分別為自然邊界條件%插值點(diǎn)x的坐標(biāo)：x0symst;f=0.0;f0=0.0;if(length(x)==length(y))n=length(x);elsedisp('x和y的維數(shù)不相等');return;endfori=1:nif(x(i)<=x0)&&(x(i+1)>=x0)index=i;break;endendA=diag(2*ones(1,n));A(1,2)=1;A(n,n-1)=1;u=zeros(n-2,1);lamda=zeros(n-1,1);c=zeros(n,1);fori=2:n-1u(i-1)=(x(i)-x(i-1))/(x(i+1)-x(i-1));lamda(i)=(x(i+1)-x(i))/(x(i+1)-x(i-1));c(i)=3*lamda(i)*(y(i)-y(i-1))/(x(i)-x(i-1))+3*u(i-1)*(y(i+1)-y(i))/(x(i+1)-x(i));A(i,i+1)=u(i-1);A(i,i-1)=lamda(i);endc(1)=3*(y(2)-y(1))/(x(2)-x(1))-(x(2)-x(1))*y2_1/2;c(n)=3*(y(n)-y(n-1))/(x(n)-x(n-1))-(x(n)-x(n-1))*y2_N/2;m=zgf(A,c);h=x(index+1)-x(index);f=y(index)*(2*(t-x(index))+h)*(t-x(index+1))^2/h/h/h+y(index+1)*(2*(x(index+1)-t)+h)*(t-x(index))^2/h/h/h+m(index)*(t-x(index))*(x(index+1)-t)^2/h/h-m(index+1)*(x(index+1)-t)*(t-x(index))^2/h/h;f0=subs(f,'t',x0);其中的zgf函數(shù)為追趕法，其程序?yàn)閒unctionx=zgf(A,b)n=rank(A);for(i=1:n)if(A(i,i)==0)disp('Error:對(duì)角有元素為0！');return;endend;d=ones(n,1);a=ones(n-1,1);c=ones(n-1);for(i=1:n-1)a(i,1)=A(i+1,i);c(i,1)=A(i,i+1);d(i,1)=A(i,i);endd(n,1)=A(n,n);for(i=2:n)d(i,1)=d(i,1)-(a(i-1,1)/d(i-1,1))*c(i-1,1);b(i,1)=b(i,1)-(a(i-1,1)/d(i-1,1))*b(i-1,1);endx(n,1)=b(n,1)/d(n,1);for(i=(n-1):-1:1)x(i,1)=(b(i,1)-c(i,1)*x(i+1,1))/d(i,1);end在MATLAB中輸入指令>>[f,f0]=scyt(x,y,0,0)f=f0=得三次樣條插值函數(shù)S(x)=1000/729*(27/5*x-1377/100)*(x-39/10)^2+1000/729*(522/25-24/5*x)*(xxxx)*(x-3)^2>>xi=0.9:0.01:13.3;yi=interp1(x,y,xi,'spline');>>title('試驗(yàn)一--三次樣條插值圖示')〔2〕用拉格朗日法插值%定義Lagrange程序functionf=Language(x,y,x0)symst;if(length(x)==length(y))n=length(x);elsedisp('x和y的維數(shù)不相等');return;endf=0.0;for(i=1:n)l=y(i);for(j=1:i-1)l=l*(t-x(j))/(x(i)-x(j));end;for(j=i+1:n)l=l*(t-x(j))/(x(i)-x(j));end;f=f+l;simplify(f);if(i==n)if(nargin==3)f=subs(f,'t',x0);elsef=collect(f);f=vpa(f,6);endendend>>Language(x,y)ans=52462.6*t+189995.*t^3-189851.*t^4+136778.*t^5-11.3161*t^12-.277283e-6*t^18+1.18284*t^13-73866.6*t^6+.111076e-4*t^17-.976904e-1*t^14+.427949e-8*t^19-.307453e-10*t^20+30677.6*t^7+2564.20*t^9-9968.98*t^8+.628590e-2*t^15-525.813*t^10-9652.78-.308159e-3*t^16+86.2514*t^11-128683.*t^2二、Wilson矩陣，且向量，那么方程組有準(zhǔn)確解。⑴用Matlab內(nèi)部函數(shù)求，的所有特征值和；⑵令，解方程組，并求出向量和，從理論結(jié)果和實(shí)際計(jì)算結(jié)果兩方面分析方程組解的相對(duì)誤差與的相對(duì)誤差的關(guān)系；⑶再改變擾動(dòng)矩陣〔其元素的絕對(duì)值不超過0.005〕，重復(fù)第2問。解：〔1〕A=[10787;7565;86109;75910];b=[32;23;33;31];M=det(A)M=1A的所以特征值：>>D=eig(A)D=0.01020.84313.858130.2887條件數(shù)>>n=30.2887/0.0102n=2.9695e+003所以A的行列式為1，cond(A)2=2.9695e+003>>B=[1078.17.2;7.085.0465;85.989.899;6.994.9999.98];>>b=[32;23;33;31];>>[rank(B),rank([B,b])]ans=44>>x=B\bx=-5.476111.4934-1.42922.4838求假設(shè)X=>>x=B\b;x1=[1;1;1;1];X=x-x1X=-6.476110.4934-2.42921.4838求>>norm(X)ans=12.655212.655就是。%求>>norm(X)/norm(x1)ans=6.32766.3276即為>>norm(a)ans=0.2244>>norm(A)ans=30.2887>>norm(a)/norm(A)%求ans=0.0074所以=0.0074inv(A)%求A的逆矩陣，下求>>d=inv(A);>>norm(d)*norm(a)*norm(x)ans=288.4858>>1-norm(d*(a))ans=-12.3693>>288.4858/-12.3693ans=-23.3227所以=-23.322>>norm(X)ans=0.0074所以>(3)改變，取a2=[00.0020.0010.003;0.0010.00400.001;0.003-0.004-0.0010;-0.001-0.0020-0.003]B2=a2+A;C2=[Bb]b=[32;23;33;31]>>rank(B2)ans=4>>rank(C2)ans=4rank(B2)=rank(C2)所以擾動(dòng)矩陣有唯一解>>x2=B2\bx2=1.06490.88931.02720.9859>>x=B\b;x1=[1;1;1;1];X=x-x1%求〔設(shè)X=〕X2=-6.476110.4934-2.42921.4838norm(X2)%求>>norm(X2)ans=12.655212.6552就是>>norm(X)/norm(x1)%求>>norm(X2)/norm(x1)ans=6.3276所以=6.3276>>norm(a2)ans=0.0071>>norm(a)/norm(A)%求>>norm(a2)/norm(A)ans=2.3336e-004所以=2.3336e-004norm(d)*norm(a2)*norm(x2)ans=9.0875>1-norm(d*(a2))ans=0.6943norm(X2)ans=12.65529.0875/0.6943ans=13.0887所以=13.0887<三、解三對(duì)角線性方程組的追趕法及其應(yīng)用⑴編寫解三對(duì)角線性方程組的追趕法的通用程序，并應(yīng)用于方程組，檢驗(yàn)程序的正確性；〔解為〕⑵求微分方程邊值問題的數(shù)值解〔取步長〕，并與精確解比擬〔精確解為〕。說明：離散化微分方程時(shí)，解：clearall;a=[2,2,2,2,2];b=[-1,-1,-1,-1];c=[-1,-1,-1,-1];r=[1,0,0,0,0];n=length(a);b=[0,b];u(1)=r(1)/a(1);v(1)=c(1)/a(1);fork=2:n-1u(k)=(r(k)-u(k-1)*b(k))/(a(k)-v(k-1)*b(k));v(k)=c(k)/(a(k)-v(k-1)*b(k));endu(n)=(r(n)-u(n-1)*b(n))/(a(n)-v(n-1)*b(n));x(n)=u(n);fork=n-1:-1:1x(k)=u(k)-v(k)*x(k+1);endfprintf('èy????·?3ì×éμ??a?a\n')fork=1:nfprintf('x(%1d)=%10.8f\n',k,x(k))end>>zhuiganfa%調(diào)用追趕法三對(duì)角方程組的解為x(1)=0.83333333x(2)=0.66666667x(3)=0.50000000x(4)=0.33333333x(5)=0.16666667和結(jié)果很好地吻合。公元1225年，比薩的數(shù)學(xué)家Leonardo研究了方程，得到一個(gè)根，沒有人知道他用什么方法得到這個(gè)值。對(duì)于這個(gè)方程，分別用以下方法：⑴迭代法；⑵迭代法；⑶對(duì)⑴的Steffensen加速方法；⑷對(duì)⑵的Steffensen加速方法；⑸Newton法。求方程的根〔可取〕，計(jì)算到Leonardo所得到的準(zhǔn)確度。解：由題意編寫m文件如下function[x0,k,er,x]=diedai(g,x0,wucha,max)%g是給定的迭代函數(shù)%x0是給定的初始值%wucha是規(guī)定的誤差范圍%max是所應(yīng)許的最大迭代次數(shù)%k是迭代次數(shù)加1%x是不動(dòng)點(diǎn)近似值%x〔x1，x2...，xn〕X(1)=x0;fork=2:maxX(k)=feval('g',X(k-1));k,er=abs(X(k)-X(k-1))x=X(k);if(er<wucha),break;endifk==maxdisp('超出迭代次數(shù)');endend其中定義的g函數(shù)是functiony=g(x);y=20/(x^2+2*x+10);在命令窗口中輸入>>diedai('g',1,10^(-9),15)k=15er=超出迭代次數(shù)X=Columns1through31.000000000000001.538461538461541.29501915708812Columns4through6Columns7through91.365929788170651.370086003401821.36824102361284Columns10through121.369059812007481.368696397555521.36885768862873Columns13through151.368786102577991.368817874396091.36880377314363ans=1取解為1.36880377314363對(duì)于，只需修改diedai.m文件中的g，把其改為g1，編寫m文件functiony=g1(x);y=(20-2*x^2-x^3)/10;在命令窗口中輸入diedai('g1',1,10^(-9),30)…..k=24er=1.36601568861328k=25er=1.36860974051282k=26er=1.36942327571766取解為1.36860974051282牛頓法：編寫m文件function[p1,er,k,y]=ndf(f,df,p0,tol,max)%f是非線性函數(shù)%df是f的微商%p0是初始值%tol是給定的允許誤差%max是迭代的最大次數(shù)%p1是牛頓法求得的近似解%er是p0的誤差估計(jì)%k是迭代次數(shù)%y=f〔p1〕p0,feval('f',p0)fork=1:maxp1=p0-feval('f',p0)/feval('df',p0);er=abs(p1-p0);p0=p1;p1,er,k,y=feval('f',p1)if((er<tol)||(y==0)),break,endend定義函數(shù)m文件functiony=f(x)y=x^3+2*x^2+10*x-20;定義微商函數(shù)m文件functiony=df(x)y=3*x^2+4*x+10;在命令窗口輸入>>ndf('f','df',1,10^(-9),10)p0=k=11ans=-7p1=1.36933647058824er=er=0.04242823529412k=k=24y=y=3.907985046680551e-014p1=p1=1.368808188617531.36880810782137K=er=31.776356839400251e-015y=k=1.704487072373695e-0065p1=y=1.368808107821370er=ans=1.36880810782137由結(jié)果知道牛頓法迭代到第三次已經(jīng)到達(dá)所要求的精度故方程的根為1.368808107821373.編寫Steffensen.m文件i=2;x0=1;%設(shè)初始值f=inline('20/(x^2+2*x+10)');%迭代格式y(tǒng)=f(x0);z=f(y);x1=x0-(y-x0).^2/(z-2*y+x0);S.result=[x0;x1];whileabs(x1-x0)>=1e-9%迭代精度x0=x1;y=f(x0);z=f(y);x1=x0-(y-x0).^2/(z-2*y+x0);i=i+1;S.result(i)=x1;endS.step=[(0:i-1)]';fprintf('迭代步數(shù)為:\t%d\n',i-1);forj=1:ifprintf('%10d',S.step(j));fprintf('\t');f在命令窗口輸入Steffensen迭代步數(shù)為: 40 1.00000001 1.37081392 1.36880823 1.36880814 1.3688081分析結(jié)果知，Steffensen迭代加速步數(shù)減少了，到第四步已經(jīng)到達(dá)了精度要求。4.把迭代格式改為(20-2*x^2-x^3)/10，保存，在命令窗口輸入Steffensen迭代步數(shù)為: 50 1.00000001 1.33349212 1.36841543 1.36880814 1.36880815 1.3688081分析結(jié)果知，Steffensen迭代加速步數(shù)減少了，到第五步已經(jīng)到達(dá)了精度要求。五、用不同的數(shù)值方法計(jì)算積分的近似值，其中⑴取不同的步長，分別用復(fù)合梯形公式和復(fù)合辛普森公式計(jì)算積分，比擬兩個(gè)公式的計(jì)算效果，是否存在一個(gè)最小的，使得精度不能再被改善？⑵用龍貝格求積公式，取，并打印出T-表。解：(1)用復(fù)合梯形公式，編寫fhtx.m文件functions=fhtx(f,a,b,n)%f是被積函數(shù)%ab分別為積分的上下限%n是子區(qū)間的個(gè)數(shù)%s是梯形總面積h=(b-a)/n;s=0;fork=1:(n-1)x=a+k*h;s=s+feval('f',x);ends=h*(feval('f',a)+feval('f',b))/2+h*s;編寫被積函數(shù)文件tf.mfunctiony=f(x);y=(10/x)^2*sin(10/x);在命令窗口輸入>>fhtx('tf',1,3,10)ans=-4.7789>>fhtx('tf',1,3,15)ans=-2.8604>>fhtx('tf',1,3,20)ans=-2.2208用復(fù)化辛普森公式，編寫fhxps.m文件functions=fhxps(f,a,b,n)%f是被積函數(shù)%ab分別為積分的上下限%n是子區(qū)間的個(gè)數(shù)%s是梯形總面積h=(b-a)/(2*n);s1=0;s2=0;fork=1:nx=a+(2*k-1)*h;s1=s1+feval('f',x);endfork=1:n-1x=a+2*k*h;s2=s2+feval('f',x);ends=h*(feval('f',a)+feval('f',b)+4*s1+2*s2)/3;在命令窗口輸入>>fhxps('f',1,3,15)ans=-1.4136>>fhxps('f',1,3,20)ans=-1.4220取n較大時(shí)>>fhxps('f',1,3,1000)ans=-1.42602475569236>>fhtx('tf',1,3,1000)ans=-1.42633620719670>>fhtx('tf',1,3,2000)ans=>>fhxps('f',1,3,2000)ans-1.42