MATLAB上机题

1、信计、数应专业10级数值方法计算实习题要求：1、用Matlab语言或你熟悉的其他语言编写程序，使之尽可能具有通用性；2、根据上机计算实践，对所使用的数值方法的特点、性质、有效性、误差和收敛性等方面进行必要的讨论和分析；3、完成计算后写出实验报告，内容包括：课题名称、解决的问题、采用的数值方法、算法程序、数值结果、对实验结果的讨论和分析等；4、特别说明：严禁抄袭，否则一经发现，所有雷同实验报告最多评为及格。1、 在区间-1,1上分别取用两组等距节点对龙格函数作多项式插值及三次样条插值，对每个值，分别画出插值函数即的图形。解：（1）多项式插值先建立一个多项式插值的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当n=10时，我们在命令窗口中输入以下的命令：clear,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')得到插值函数和f（x）图形：当n=20时，我们在命令窗口中输入以下的命令：clear,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')得到插值函数和f（x）图形：（2）三次样条插值先建立一个多项式插值的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当n=10时，我们在命令窗口中输入以下的命令：clear,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,'.')结果如图：当n=20时，我们在命令窗口中输入以下的命令：clear,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,'.')结果如图：二、已知Wilson矩阵，且向量，则方程组有准确解。用Matlab 内部函数求，的所有特征值和；令，解方程组，并求出向量和，从理论结果和实际计算结果两方面分析方程组解的相对误差与的相对误差的关系；再改变扰动矩阵（其元素的绝对值不超过0.005），重复第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-范数条件数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求 假设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)所以扰动矩阵有唯一解>> x2=B2bx2 = 1.0649 0.8893 1.02720.9859>> x=Bb;x1=1;1;1;1;X=x-x1 %求（设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， <三、解三对角线性方程组的追赶法及其应用编写解三对角线性方程组的追赶法的通用程序，并应用于方程组，检验程序的正确性；（解为）求微分方程边值问题的数值解（取步长），并与精确解比较（精确解为）。说明：离散化微分方程时，解：（1）编写解三对角线性方程组的追赶法的通用程序：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('三对角方程组的解为：n')for k=1:n fprintf('x(%1d)=%10.8fn',k,x(k)end>> zhuiganfa %调用追赶法三对角方程组的解为：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('解曲线','解点','粗略解')指令窗口:sinit = x: 0 1 2 3 y: 2x4 doublesol = x: 1x17 double y: 2x17 double yp: 2x17 double solver: 'bvp4c'四、公元1225年，比萨的数学家Leonardo研究了方程，得到一个根，没有人知道他用什

18、么方法得到这个值。对于这个方程，分别用下列方法：迭代法；迭代法；对的Steffensen加速方法；对的Steffensen加速方法；Newton法。求方程的根（可取），计算到Leonardo所得到的准确度。解：(1)由题意编写m文件如下：functionx0,k,er,x=diedai(g,x0,wucha,max)%g是给定的迭代函数%x0是给定的初始值%wucha是规定的误差范围%max是所应许的最大迭代次数%k是迭代次数加1%x是不动点近似值%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('超出迭代次数'); end end %其中定义的g函数是：function y=g(x); y=20/(x2+2*x+10);在命令窗口中输入:>> diedai('g',1,10(-9),15)k = 15er = 1.410125245193683e-005超出迭代次数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）对于 ，只需修改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;%设初始值f=inline('20/(x2+2*x+10)');%迭代格式y=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('迭代步数为:t%dn',i-1);for j=1:i fprintf('%10d',S.step(j);fprintf('t'); f在命令窗口输入Steffensen迭代步数为:4 0 1.0000000 1 1.3708139 2 1

23、.3688082 3 1.3688081 4 1.3688081分析结果知，Steffensen迭代加速步数减少了，到第四步已经达到了精度要求。（4）把迭代格式改为(20-2*x2-x3)/10，保存，在命令窗口输入Steffensen迭代步数为:5 0 1.0000000 1 1.3334921 2 1.3684154 3 1.3688081 4 1.3688081 5 1.3688081分析结果知Steffensen迭代加速步数减少了，到第五步已经达到了精度要求。（5）牛顿法：编写m文件functionp1,er,k,y=ndf(f,df,p0,tol,max)%f是非线性函数%df是f的

24、微商%p0是初始值%tol是给定的允许误差%max是迭代的最大次数%p1是牛顿法求得的近似解%er是p0的误差估计%k是迭代次数%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%定义函数m文件function y=f(x)y=x3+2*x2+10*x-20;%定义微商函数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由结果知道牛顿法迭代到第三次已经达到所要求的精度故方程的根为1.368808107821373五、用不同的数值方法计算积分的近似值，其中 取不同的步长，分别用复合梯形公式和复合辛普森公式计算积分，比较两个公式的计算效果，是否存在一个最小

27、的，使得精度不能再被改善？ 用龙贝格求积公式，取，并打印出T-表。解：（1）复合梯形公式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复合辛普森公式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