MATLAB上機(jī)題

1、信計(jì)、數(shù)應(yīng)專業(yè)10級(jí)數(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)為及格。1、 在區(qū)間-1,1上分別取用兩組等距節(jié)點(diǎn)對(duì)龍格函數(shù)作多項(xiàng)式插值及三次樣條插值，對(duì)每個(gè)值，分別畫出插值函數(shù)即的圖形。解：（1）多項(xiàng)式插值先建立一個(gè)多項(xiàng)式插值的M-file；輸入如下的命令（如牛

2、頓插值公式）：function C,D=newpoly(X,Y)n=length(X);D=zeros(n,n) D(:,1)=Y' for j=2:n for k=j:n D(k,j)=(D(k,j-1)- D(k-1,j-1)/(X(k)-X(k-j+1); endendC=D(n,n);for k=(n-1):-1:1 C=conv(C,poly(X(k) m=length(C); C(m)= C(m)+D(k,k);end當(dāng)n=10時(shí)，我們?cè)诿畲翱谥休斎胍韵碌拿睿篶lear,clf,hold on; X=-1:0.2:1; Y=1./(1+25*X.2); C,D=newp

3、oly(X,Y); x=-1:0.01:1; y=polyval(C,x); plot(x,y,X,Y,'.'); grid on; xp=-1:0.2:1; z=1./(1+25*xp.2); plot(xp,z,'r')得到插值函數(shù)和f（x）圖形：當(dāng)n=20時(shí)，我們?cè)诿畲翱谥休斎胍韵碌拿睿篶lear,clf,hold on; X=-1:0.1:1; Y=1./(1+25*X.2); C,D=newpoly(X,Y); x=-1:0.01:1; y=polyval(C,x); plot(x,y,X,Y,'.'); grid on; xp=-

4、1:0.1:1; z=1./(1+25*xp.2); plot(xp,z,'r')得到插值函數(shù)和f（x）圖形：（2）三次樣條插值先建立一個(gè)多項(xiàng)式插值的M-file；輸入如下的命令:function S=csfit(X,Y,dx0,dxn) N=length(X)-1;H=diff(X); D=diff(Y)./H;A=H(2:N-1);B=2*(H(1:N-1)+H(2:N); C=H(2:N); U=6*diff(D);B(1)=B(1)-H(1)/2;U(1)=U(1)-3*(D(1);B(N-1)=B(N-1)-H(N)/2;U(N-1)=U(N-1)-3*(-D(N);

5、 for k=2:N-1 temp=A(k-1)/B(k-1); B(k)=B(k)-temp*C(k-1); U(k)=U(k)-temp*U(k-1); end M(N)=U(N-1)/B(N-1);for k=N-2:-1:1 M(k+1)=(U(k)-C(k)*M(k+2)/B(k);endM(1)=3*(D(1)-dx0)/H(1)-M(2)/2;M(N+1)=3*(dxn-D(N)/H(N)-M(N)/2;for k=0:N-1 S(k+1,1)=(M(k+2)-M(k+1)/(6*H(k+1); S(k+1,2)=M(k+1)/2; S(k+1,3)=D(k+1)-H(k+1)*

6、(2*M(k+1)+M(k+2)/6; S(k+1,4)=Y(k+1);end當(dāng)n=10時(shí)，我們?cè)诿畲翱谥休斎胍韵碌拿睿篶lear,clcX=-1:0.2:1;Y=1./(25*X.2+1); dx0= 0.0739644970414201;dxn= -0.0739644970414201;S=csfit(X,Y,dx0,dxn) x1=-1:0.01:-0.5;y1=polyval(S(1,:),x1-X(1);x2=-0.5:0.01:0;y2=polyval(S(2,:),x2-X(2); x3=0:0.01:0.5; y3=polyval(S(3,:),x3-X(3);x4=0.5

7、:0.01:1;y4=polyval(S(4,:),x4-X(4); plot(x1,y1,x2,y2,x3,y3,x4,y4, X,Y,'.')結(jié)果如圖：當(dāng)n=20時(shí)，我們?cè)诿畲翱谥休斎胍韵碌拿睿篶lear,clcX=-1:0.1:1;Y=1./(25*X.2+1); dx0= 0.0739644970414201;dxn= -0.0739644970414201;S=csfit(X,Y,dx0,dxn) x1=-1:0.01:-0.5;y1=polyval(S(1,:),x1-X(1);x2=-0.5:0.01:0;y2=polyval(S(2,:),x2-X(2);

8、x3=0:0.01:0.5; y3=polyval(S(3,:),x3-X(3);x4=0.5:0.01:1;y4=polyval(S(4,:),x4-X(4); plot(x1,y1,x2,y2,x3,y3,x4,y4, X,Y,'.')結(jié)果如圖：二、已知Wilson矩陣，且向量，則方程組有準(zhǔn)確解。用Matlab 內(nèi)部函數(shù)求，的所有特征值和；令，解方程組，并求出向量和，從理論結(jié)果和實(shí)際計(jì)算結(jié)果兩方面分析方程組解的相對(duì)誤差與的相對(duì)誤差的關(guān)系；再改變擾動(dòng)矩陣（其元素的絕對(duì)值不超過0.005），重復(fù)第2問。解：（1）>> A=10 7 8 7;7 5 6 5;8 6 1

9、0 9;7 5 9 10;b=32;23;33;31;>> C=det(A) %求A的行列式C = 1>> X=eig(A) %求A的所有特征值X = 0.0102 0.8431 3.8581 30.2887>> D=cond(A) %求A的2-范數(shù)條件數(shù)D = 2.9841e+003所以A的行列式為1，cond(A)2=2.9695e+003(2)>> B=10 7 8.1 7.2;7.08 5.04 6 5;8 5.98 9.89 9;6.99 4.99 9 9.98;>> b=32;23;33;31;>> rank(

10、B),rank(B,b)ans = 4 4>> x=Bb x = -5.4761 11.4934 -1.4292 2.4838求 假設(shè)X= >> x=Bb;x1=1;1;1;1;X=x-x1 X = -6.4761 10.4934 -2.42921.4838%求 >>norm(X) ans = 12.655212.655就是。 %求 >> norm(X)/norm(x1)ans = 6.3276 6.3276 即為>> norm(a)ans = 0.2244>> norm(A)ans = 30.2887>> n

11、orm(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=0 0.002 0.001 0.003;0.001 0.004 0 0.001;0.003 -0.00

12、4 -0.001 0;-0.001 -0.002 0 -0.003B2=a2+A; C2=B b b=32;23;33;31>> rank(B2) ans = 4>> rank(C2)ans = 4rank(B2)= rank(C2)所以擾動(dòng)矩陣有唯一解>> x2=B2bx2 = 1.0649 0.8893 1.02720.9859>> x=Bb;x1=1;1;1;1;X=x-x1 %求（設(shè)X= ）X2 = -6.4761 10.4934 -2.4292 1.4838norm(X2) %求 >> norm(X2)ans = 12.65

13、5212.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-004>> norm(d)*norm(a2)*norm(x2)ans = 9.0875>> 1-norm(d*(a2)ans = 0.6943>>norm(X2)ans =

14、 12.6552 9.0875/0.6943ans = 13.0887所以= 13.0887， <三、解三對(duì)角線性方程組的追趕法及其應(yīng)用編寫解三對(duì)角線性方程組的追趕法的通用程序，并應(yīng)用于方程組，檢驗(yàn)程序的正確性；（解為）求微分方程邊值問題的數(shù)值解（取步長），并與精確解比較（精確解為）。說明：離散化微分方程時(shí)，解：（1）編寫解三對(duì)角線性方程組的追趕法的通用程序：clear all;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);for k=

15、2:n-1 u(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);for k=n-1:-1:1 x(k)=u(k)-v(k)*x(k+1);endfprintf('三對(duì)角方程組的解為：n')for k=1:n fprintf('x(%1d)=%10.8fn',k,x(k)end>> zhuiganfa %調(diào)用追趕法三對(duì)角方程組的解為：x(1)=0.833333

16、33x(2)=0.66666667x(3)=0.50000000x(4)=0.33333333x(5)=0.16666667（2）邊值解M-文件：sinit=bvpinit(0:pi,2;0)odefun=inline('-y(2);abs(y(2)-y(1)+exp(t)-3*sin(t)','t','y');bcfun=inline('ya(1)+2;yb(1)-exp(pi)-3','ya','yb');sol=bvp4c(odefun,bcfun,sinit)t=linspace(0,pi,

17、101);y=deval(sol,t);plot(t,y(2,:),sol.x,sol.y(2,:),'o',sinit.x,sinit.y(2,:),'s')legend('解曲線','解點(diǎn)','粗略解')指令窗口:sinit = x: 0 1 2 3 y: 2x4 doublesol = x: 1x17 double y: 2x17 double yp: 2x17 double solver: 'bvp4c'四、公元1225年，比薩的數(shù)學(xué)家Leonardo研究了方程，得到一個(gè)根，沒有人知道他用什

18、么方法得到這個(gè)值。對(duì)于這個(gè)方程，分別用下列方法：迭代法；迭代法；對(duì)的Steffensen加速方法；對(duì)的Steffensen加速方法；Newton法。求方程的根（可?。?jì)算到Leonardo所得到的準(zhǔn)確度。解：(1)由題意編寫m文件如下：functionx0,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;for k=2:max X(k)=feval('g',X(k-1); k,er=abs(X

19、(k)-X(k-1)x=X(k);if(er<wucha), break;endif k=max disp('超出迭代次數(shù)'); end end %其中定義的g函數(shù)是：function y=g(x); y=20/(x2+2*x+10);在命令窗口中輸入:>> diedai('g',1,10(-9),15)k = 15er = 1.410125245193683e-005超出迭代次數(shù)X =Columns 1 through 4 1.00000000000000 1.53846153846154 1.29501915708812 1.4018253

20、0944860Columns 5 through 8 1.35420939040429 1.37529809248738 1.36592978817065 1.37008600340182Columns 9 through 12 1.36824102361284 1.36905981200748 1.36869639755552 1.36885768862873Columns 13 through 15 1.36878610257799 1.36881787439609 1.36880377314363ans =1取解為1.36880377314363（2）對(duì)于 ，只需修改diedai.m文件

21、中的g，把其改為g1，編寫m文件function y=g1(x);y=(20-2*x2-x3)/10;在命令窗口中輸入diedai('g1',1,10(-9),30).k = 24er = 1.36601568861328k = 25er = 1.36860974051282k = 26er = 1.36942327571766取解為1.36860974051282（3）編寫Steffensen.m文件i=2;x0=1;%設(shè)初始值f=inline('20/(x2+2*x+10)');%迭代格式y(tǒng)=f(x0);z=f(y);x1=x0-(y-x0).2/(z-2*

22、y+x0);S.result=x0;x1;while abs(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%dn',i-1);for j=1:i fprintf('%10d',S.step(j);fprintf('t'); f在命令窗口輸入Steffensen迭代步數(shù)為:4 0 1.0000000 1 1.3708139 2 1

23、.3688082 3 1.3688081 4 1.3688081分析結(jié)果知，Steffensen迭代加速步數(shù)減少了，到第四步已經(jīng)達(dá)到了精度要求。（4）把迭代格式改為(20-2*x2-x3)/10，保存，在命令窗口輸入Steffensen迭代步數(shù)為:5 0 1.0000000 1 1.3334921 2 1.3684154 3 1.3688081 4 1.3688081 5 1.3688081分析結(jié)果知Steffensen迭代加速步數(shù)減少了，到第五步已經(jīng)達(dá)到了精度要求。（5）牛頓法：編寫m文件functionp1,er,k,y=ndf(f,df,p0,tol,max)%f是非線性函數(shù)%df是f的

24、微商%p0是初始值%tol是給定的允許誤差%max是迭代的最大次數(shù)%p1是牛頓法求得的近似解%er是p0的誤差估計(jì)%k是迭代次數(shù)%y=f（p1）p0,feval('f',p0)for k=1:max p1=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文件function y=f(x)y=x3+2*x2+10*x-20;%定義微商函數(shù)m文件fu

25、nction y=df(x)y=3*x2+4*x+10;%在命令窗口輸入>> ndf('f','df',1,10(-9),10)p0 = 1ans = -7p1 = 1.411764705882353er = 0.411764705882353k = 1y = 0.917565642173823p1 = 1.369336470588235er = 0.042428235294118k = 2y = 0.011148119412447p1 = 1.368808188617532er = 5.282819707033148e-004k = 3y = 1.

26、704487072373695e-006p1 = 1.368808107821375er = 8.079615732015100e-008k = 4y = 3.907985046680551e-014p1 = 1.368808107821373er = 1.776356839400251e-015k = 5y = 0ans = 1.368808107821373由結(jié)果知道牛頓法迭代到第三次已經(jīng)達(dá)到所要求的精度故方程的根為1.368808107821373五、用不同的數(shù)值方法計(jì)算積分的近似值，其中 取不同的步長，分別用復(fù)合梯形公式和復(fù)合辛普森公式計(jì)算積分，比較兩個(gè)公式的計(jì)算效果，是否存在一個(gè)最小

27、的，使得精度不能再被改善？ 用龍貝格求積公式，取，并打印出T-表。解：（1）復(fù)合梯形公式M文件：function I,step = CombineTraprl(f,a,b,eps)if(nargin =3) eps=1.0e-4;endn=1;h=(b-a)/2;I1=0;I2=(subs(sym(f),findsym(sym(f),a)+subs(sym(f),findsym(sym(f),b)/h;while abs(I2-I1)>eps n=n+1; h=(b-a)/n; I1=I2; I2=0; for i=0:n-1 x=a+h*i; x1=x+h; I2=I2+(h/2)*(

28、subs(sym(f),findsym(sym(f),x)+subs(sym(f),findsym(sym(f),x1); endendI=I2;step=n;在命令窗口輸入： >> q s= CombineTraprl('(10/x)2*sin(10/x)',1,3,1.0e-3)q = -1.4682s =86復(fù)合辛普森公式M文件：function I,step = IntSimpson(f,a,b,type,eps)if(type=3 && nargin=4) eps=1.0e-4; endI=0;switch type case 1, I=(

29、b-a)/6)*(subs(sym(f),findsym(sym(f),a)+. 4*subs(sym(f),findsym(sym(f),(a+b)/2)+. subs(sym(f),findsym(sym(f),b); step=1; case 2, I=(b-a)/8)*(subs(sym(f),findsym(sym(f),a)+. 3*subs(sym(f),findsym(sym(f),(2*a+b)/3)+ . 3*subs(sym(f),findsym(sym(f),(a+2*b)/3)+subs(sym(f),findsym(sym(f),b); step=1; case 3, n=2; h=(b-a)/2; I1=0; I2=(subs(sym(f),findsym(sym(f),a)+subs(sym(f),findsym(sym(f),b)/h; while abs(I2-I1)>eps n=n+1; h=(b-a)/n; I1=I2; I2=0; for i=0:n-1 x=a+h*i; x1=x+h; I2=I2+(h/6)*(subs(sym(f),findsym(sym(f),x)+. 4*subs(sym(f),findsym(sym(f),(x+x1)/2)+. subs(sym(f),findsym(sym(f),x1); end