数值分析 (2)

1、1.1程序：n=input('输入n '); for x=1:n p=5*x+1; for k=1:5 p=5*p/4+1; end if p=fix(p) break endenddisp(x,p)结果：输入n：3000 1023 156211.3p=0 0;10 0; %p为初始两个点的坐标，第一列为x坐标，二列为yn=2; %n为节点数A=cos(pi/3) -sin(pi/3);sin(pi/3) cos(pi/3); %旋转矩阵for k=1:10 d=diff(p)/3; %两相临点坐标差，得到相邻两点的向量 %d为向量的1/3与题中三等分对应 m=4*n-3; %

2、迭代式 q=p(1:n-1,:); %原点与前n-1个点形成向量 p(5:4:m,:)=p(2:n,:); %迭代后处于4k+1位置上的点的坐标为迭代前的相应坐标p(2:4:m,:)=q+d; %用向量方法计算迭代后处于4k+2位置上的点的坐标 p(3:4:m,:)=q+d+d*A' %用向量方法计算迭代后处于4k+3位置上的点的坐标 p(4:4:m,:)=q+2*d; %迭代后新的结点数目n=m; %迭代后新的结点数目end plot(p(:,1),p(:,2) %绘出每相邻两个点的连线axis(0 10 0 10)2.1（1）将给定的点带入得方程组：整理后=>AX=R（2）在

3、命令窗口输入：A=536052 53605*12052 60262 107210 12052; 58460 58460*22358 111792 116920 22358; 628592 62859*33908 169542 125718 33908; 666622 66662*46984 234922 133334 46984; 688942 688942*2 688942 68894*2 68894*2;b=-1;-1;-1;-1;-1;Ab结果：ans = 1.0e-004 * -0.000000000024637 -0.000000848114996 -0.00000241202170

4、42.3（1）在matlab命令窗口输入程序： A=20 2 3;1 8 1;2 -3 15; L U=lu(A)输出结果：L = 1.0000 0 0 0.0500 1.0000 0 0.1000 -0.4051 1.0000U = 20.0000 2.0000 3.0000 0 7.9000 0.8500 0 0 15.0443(2) 在matlab命令窗口输入程序：iiA=20 2 3;1 8 1;2 -3 15;b1=1 0 0'b2=0 1 0'b3=0 0 1' L U=lu(A);结果：X1=U(Lb1)X1 = 0.0517 -0.0055 -0.008

5、0 X2=U(Lb2)X2 = -0.0164 0.1237 0.0269 X3=U(Lb3)X3 = -0.0093 -0.00720.0665输入：inv(A)ans = 0.0517 -0.0164 -0.0093 -0.0055 0.1237 -0.0072 -0.0080 0.0269 0.0665(3) 在matlab命令窗口输入程序： A=4 2 -3 -1 2 1 0 0 0 0;8 6 -5 -3 6 5 0 1 0 0;4 2 -2 -1 3 2 -1 0 3 1;0 -2 1 5 -1 3 -1 1 9 4;-4 2 6 -1 6 7 -3 3 2 3;8 6 -8 5

6、7 17 2 6 -3 5;0 2 -1 3 -4 2 5 3 0 1;16 10 -11 -9 17 34 2 -1 2 2;4 6 2 -7 13 9 2 0 12 4;0 0 -1 8 -3 -24 -8 6 3 -1; B=5 12 3 2 3 46 13 38 19 -21'L U=lu(A);X=U(LB)结果：X = 1.0000 -1.0000 0.0000 1.0000 2.0000 0.0000 3.0000 1.0000 -1.0000 2.00003.1Sinx的taylor展开式的前三项y,y1,y2输入程序：x=0:0.01:pi;y0=sin(x);y=x

7、;y1=x-x.3/6;y2=x-x.3/6+x.5/120;subplot(2,2,1);plot(x,y0);title('sinx');grid on;subplot(2,2,2);plot(x,y,'g');title('x');grid on;subplot(2,2,3);plot(x,y1,'r');title('x-x.3/6');grid on;subplot(2,2,4);plot(x,y2,'y');title('x-x.3/6+x.5/120');grid on

8、;X为0,pi时的图像逼近：X为0,pi/2时的图像逼近：3.2程序： syms x1 x2 y1 y2z=3*x12+3*y12+3*x22+3*y22-6*y1-12*y2-36*x2-2*x1*x2-2*y1*y2+199dx1=diff(z,'x1')dx2=diff(z,'x2')dy1=diff(z,'y1')dy2=diff(z,'y2')x1,y1,x2,y2=solve('6*x1-2*x2=0','6*x2-36-2*x1=0','6*y1-6-2*y2=0',&

9、#39;6*y2-12-2*y1=0','x1','y1','x2','y2')结果：dx1 =6*x1 - 2*x2dx2 =6*x2 - 2*x1 - 36 dy1 =6*y1 - 2*y2 - 6 dy2 =6*y2 - 2*y1 - 12 x1 =9/4y1 =27/4 x2 =15/8y2 =21/84.1程序：x=-6:0.1:6;y1=sin(x);y2=1./x;plot(x,y1);hold on;plot(x,y2,'r');画图：观察图形可得:方程的实根分布在1,1.5,2.5,3两个

10、区间,接近0点时没有交点,估计在-6,0.2区间-1,-1.5,-3,-2.5之间,接近零点也没有交点,两者分散.4.3（1，） x=1 -55 1320 -18150 157773 -902055 3416930 -8409500 12753576 -10628640 6328800;y=roots(x)结果：>> x=1 -55 1320 -18150 157773 -902055 3416930 -8409500 12753576 -10628640 6328800;y=roots(x)y = 10.6051 + 1.0127i 10.6051 - 1.0127i 8.585

11、0 + 2.7898i 8.5850 - 2.7898i 5.5000 + 3.5058i 5.5000 - 3.5058i 2.4150 + 2.7898i 2.4150 - 2.7898i 0.3949 + 1.0127i 0.3949 - 1.0127i（2，）（3，）>> x=1 -56 1320 -18150 157773 -902055 3416930 -8409500 12753576 -10628640 6328800;y=roots(x)y = 21.7335 7.3501 + 7.7973i 7.3501 - 7.7973i 4.3589 + 3.3285i 4

12、.3589 - 3.3285i 5.1831 2.4378 + 2.7974i 2.4378 - 2.7974i 0.3949 + 1.0127i 0.3949 - 1.0127i根的变化挺显著的，说明多项式系数微小的变化可以引起根的明显变化！5.1(1,)lagrange 差值！function f=lag(x,y,x0)syms t;if length(x)=length(y) n=length(x);else disp('x yάÊý²»Í¬'); return;endf=0.0;for i

13、=1:n l=y(i); for j=1:i-1 l=1*(t-x(j)/(x(i)-x(j); end for j=i+1:n l=1*(t-x(j)/(x(i)-x(j); end f=f+l; simplify(f); if i=n if nargin=3 f=subs(f,'t',x0); else f=collect(f); f=vpa(f,5); end endend结果：x=0.4,0.55,0.65,0.8,0.9,1.05;y=0.41075,0.57815,0.69675,0.88811,1.02652,1.25382;f=lag(x,y,0.596)f =

14、 5.5575Newton 插值法！function f=new(x,y,x0)syms t;if length(x)=length(y) n=length(x); c(1:n)=0.0;else disp('x y άÊý²»Í¬'); return;endf=y(1);y1=0;l=1;for i=1:n-1 for j=i+1:n y1(j)=(y(j)-y(i)/(x(j)-x(i); end c(i)=y1(i+1); l=l*(t-x(i); f=f+c(i)*l; simplify(

15、f); y=y1; if i=n-1 if nargin=3 f=subs(f,'t',x0); else f=collect(f); f=vpa(f,5); end endend 结果： x=0.4,0.55,0.65,0.8,0.9,1.05;y=0.41075,0.57815,0.69675,0.88811,1.02652,1.25382;f=lag(x,y,0.596)f = 5.5575（2，）Aitken 算法！function Ait(x,y,a)n=length(x);p=ones(1,n);p(:)=y(:);for k=1:n-1 for i=k+1:n p

16、(i)=p(i)+(p(i)-p(k)/(x(i)-x(k)*(a-x(i); endendp(n)结果：x=0.4,0.55,0.65,0.8,0.9,1.05;y=0.41075,0.57815,0.69675,0.88811,1.02652,1.25382;f=Ait(x,y,0.596)f = 5.5575Neville算法！function nev(x,y,a)n=length(x);p=zeros(n);p(1,: )=y(1,: );for k=1:n-1 for i=k+1:n p(k+1,i)=(a-x(i-k)*p(k,i)-(a-x(i)*p(k,i-1)/(x(i)-x

17、(i-k); endendp(n,n)结果：x=0.4,0.55,0.65,0.8,0.9,1.05;y=0.41075,0.57815,0.69675,0.88811,1.02652,1.25382;f=lag(x,y,0.596)f = 5.5575(3,) 牛顿向前差值！function f=newf(x,y,x0)syms t;if length(x)=length(y)n=length(x);c(1:n)=0.0;else disp('x yάÊý²»Í¬'); return;endf=

18、y(1);y1=0;xx=linspace(x(1),x(n),(x(2)-x(1);if xx=x disp('²»µÈ¾à'); return;endfor i=1:n-1 for j=1:n-i y1(j)=y(j+1)-y(j); end c(i)=y1(1); l=t; for k=1:i-1 l=l*(t-k); end f=f+c(i)*l/factorial(i); simplify(f); y=y1; if i=n-1 if nargin=3 f=subs(f,'t',(x0-x(1)

19、/(x(2)-x(1); else f= colllect(f); f=vpa(f,5); end endend结果：x=0.4,0.55,0.65,0.8,0.9,1.05;y=0.41075,0.57815,0.69675,0.88811,1.02652,1.25382;f=newf(x,y,0.42)f =0.4569牛顿向后差值！function f=newb(x,y,x0)syms t;if length(x)=length(y) n=length(x); c(1:n)=0.0;else disp('x yάÊý²»

20、;ͬ'); return;endf=y(n);y1=0;xx=linspace(x(1),x(n),(x(2)-x(1);if xx=x disp('²»µÈ¾à'); return;endfor i=1:n-1 for j=i+1:n; y1(j)=y(j)-y(j-1); end c(i)=y1(n); l=t; for k=1:i-1 l=l*(t+k); end f=f+c(i)*l/factorial(i); simplify(f); y=y1; if i=n-1 if narg

21、in=3 f=subs(f,'t',(x0-x(n)/(x(2)-x(1); else f= colllect(f); f=vpa(f,5); end endend结果： x=0.4,0.55,0.65,0.8,0.9,1.05;y=0.41075,0.57815,0.69675,0.88811,1.02652,1.25382;f=newb(x,y,0.42)f = 0.56266.2function zjpfbj(f,n,a,b)for i=1:n+1 f1=eval('(x)(',f,'.*x.(i-1)'); d(i)=quad(f1,a,

22、b);endformat longH=hilb(n+1);L,U=lu(H);alpha=U(Ld');alpha结果：zjpfbj('sin(pi*x)',2,0,1)alpha = -0.050464849874539>> zjpfbj('sin(pi*x)',3,0,1)alpha = -0.050459754555635>> zjpfbj('sin(pi*x)',6,0,1)alpha = 0.254225561815206 -0.602780183279169>> zjpfbj('si

23、n(pi*x)',8,0,1)alpha = 1.0e+003 * -0.000027438509887 0.005271152811926 -0.039361877097406 0.295759697707640 2.492982607327577 -2.969493901887367 1.850869720128680 -0.469975834581631>> zjpfbj('sin(pi*x)',10,0,1)alpha = 1.0e+007 * -0.000001000675394 -0.002865115797952 0.0326354636213

24、05 0.703092988877561 -1.783426233720688作图： x=0:0.01:1;y=sin(pi*x);plot(x,y);title ('sin(pi*x)');xlabel('x');ylabel('y');x=0:0.01:1;plot(x,y);title('y1µÄͼ');xlabel('x1');ylabel('y1');x=0:0.01:1;y2plot(x,y);title('y2µÄ&

25、#205;¼');xlabel('x2');ylabel('y2'); x=0:0.01:1;plot(x,y);title('y3µÄͼ');xlabel('x3');ylabel('y3');x=0:0.01:1; plot(x,y);title('y4µÄͼ');xlabel('x4');ylabel('y4');7.3求积分问题，当u取不同的值可以得到不同的积

26、分。采用复化的辛普森公式求解。for a=0.1:0.1:3 b=4;n=10;h=(b-a)/n; p=0;q=0;f=inline('exp(-x2)/2)'); for i=0:n-1 p1=f(a+i*h+a+(i+1)*h)/2);p=p+p1;endfor j=1:n-1q1=f(a+j*h); q=q+q1; endformat longs=(1/sqrt(2*pi)*(h/6)*(f(a)+4*p+2*q+f(b)end结果： s =0.420706858216159s =0.382054652282866 s =0.344544062522814s = 0.3

27、08503267097576 s =0.096768261388249 s = 0.080724631768960s = 0.066775324333546 s = 0.054767524041190s = 0.044533769129462 s = 0.035898671181465s = 0.028684936179116 s =0.004629533589453s = 0.003435313851058 s =0.002523466314120 7.6(1)易得概率值的计算实际上是求定积分，采用复化的新普生公式进行求解：(2)实际上是求出相应的概率值满足题意要求即可。由于我们从170到194之间选出一个能满足。（3）用normrnd(170,6,10000,1)可模拟城市10