数值计算方法程序设计2.doc

《数值计算方法程序设计2.doc》由会员分享，可在线阅读，更多相关《数值计算方法程序设计2.doc（12页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。

1、,1.用matlab编写拉格朗日插值算法的程序 并且以（x=-2.00,f(x)=17.00 x=0.00,f(x)=1.00 x=1.00,f(x)=2.00 x=2.00,f（x)=17.00)为数据基础，在整个插值区间上采用拉格朗日插值算法计算f(x=0.6)，写出程序源代码，输出计算结果x0=-2.00;x1=0.00;x2=1.00;x3=2.00;y0=17.00;y1=1.00;y2=2.00;y3=17.00;x=0.6y=(x-x1).*(x-x2).*(x-x3)/(x0-x1).*(x0-x2).*(x0-x3)*y0+(x-x0).*(x-x2).*(x-x3)/(x1

2、-x0).*(x1-x2).*(x1-x3)*y1+(x-x0).*(x-x1).*(x-x3)/(x2-x0).*(x2-x1).*(x2-x3)*y2+(x-x0).*(x-x1).*(x-x2)/(x3-x0).*(x3-x1).*(x3-x2)*y3;disp(y=);disp(y);结果为：x = 0.6000y=0.25602.追赶法function x=zhuiganfa%首先说明：追赶法是适用于三对角矩阵的线性方程组求解的方法，并不适用于其他类型矩阵。%定义三对角矩阵A的各组成单元。方程为Ax=d% b为A的对角线元素(1n)，a为-1对角线元素(2n)，c为+1对角线元素(1

3、n-1)。% A=2 -1 0 0% -1 3 -2 0% 0 -2 4 -3% 0 0 -3 5a=0 -1 -2 -3;c=-1 -2 -3;b=2 3 4 5;d=6 1 -2 1;n=length(b);u0=0;y0=0;a(1)=0;%“追”的过程L(1)=b(1)-a(1)*u0;y(1)=(d(1)-y0*a(1)/L(1);u(1)=c(1)/L(1);for i=2:(n-1) L(i)=b(i)-a(i)*u(i-1); y(i)=(d(i)-y(i-1)*a(i)/L(i); u(i)=c(i)/L(i);endL(n)=b(n)-a(n)*u(n-1);y(n)=(d

4、(n)-y(n-1)*a(n)/L(n);%“赶”的过程x(n)=y(n);for i=(n-1):-1:1 x(i)=y(i)-u(i)*x(i+1);end3.特征向量的计算，幂法5.2.2 幂法的MATLAB程序用幂法计算矩阵的主特征值和对应的特征向量的MATLAB主程序function k,lambda,Vk,Wc=mifa(A,V0,jd,max1)lambda=0;k=1;Wc =1; ,jd=jd*0.1;state=1; V=V0;while(kjd)state=1;endk=k+1;Wc=Wc;endif(WcA=1 -1;2 4; V0=1,1; k,lambda,Vk,W

5、c=mifa(A,V0,0.00001,100), V,D = eig (A), Dzd=max(diag(D), wuD= abs(Dzd- lambda), wuV=V(:,2)./Vk, 运行后屏幕显示结果请注意：迭代次数k,主特征值的近似值lambda,主特征向量的近似向量Vk,相邻两次迭代的误差Wc如下：k = lambda = Wc = 33 3.00000173836804 8.691862856124999e-007Vk = V = wuV =-0.49999942054432 -0.70710678118655 0.44721359549996 -0.894428227562

6、941.00000000000000 0.70710678118655 -0.89442719099992 -0.89442719099992Dzd = wuD = 3 1.738368038406435e-006 由输出结果可看出，迭代33次，相邻两次迭代的误差Wc 8.69 19e-007，矩阵的主特征值的近似值lambda3.000 00和对应的特征向量的近似向量Vk （-0.500 00，1.000 00, lambda与例5.1.1中的最大特征值近似相等，绝对误差约为1.738 37e-006，Vk与特征向量的第1个分量的绝对误差约等于0，第2个分量的绝对值相同.由wuV可以看出，的

7、特征向量V(:,2) 与Vk的对应分量的比值近似相等.因此，用程序mifa.m计算的结果达到预先给定的精度.(2) 输入MATLAB程序 B=1 2 3;2 1 3;3 3 6; V0=1,1,1; k,lambda,Vk,Wc=mifa(B,V0,0.00001,100), V,D = eig (B), Dzd=max(diag(D), wuD= abs(Dzd- lambda), wuV=V(:,3)./Vk,运行后屏幕显示结果请注意：迭代次数k,主特征值的近似值lambda,主特征向量的近似向量Vk,相邻两次迭代的误差Wc如下：k = lambda = Wc = Dzd = wuD =

8、3 9 0 9 0Vk = wuV = 0.50000000000000 0.81649658092773 0.50000000000000 0.81649658092773 1.00000000000000 0.81649658092773V = 0.70710678118655 0.57735026918963 0.40824829046386 -0.70710678118655 0.57735026918963 0.40824829046386 0 -0.57735026918963 0.81649658092773 (3) 输入MATLAB程序 C=1 2 2;1 -1 1;4 -12

9、 1;V0=1,1,1; k,lambda,Vk,Wc=mifa(C,V0,0.00001,100), V,D = eig (C), Dzd=max(diag(D), wuD=abs(Dzd-lambda),Vzd=V(:,1),wuV=V(:,1)./Vk,运行后屏幕显示请注意：迭代次数k已经达到最大迭代次数max1,主特征值的迭代值lambda,主特征向量的迭代向量Vk,相邻两次迭代的误差Wc如下：k = lambda = Wc = 100 0.09090909090910 2.37758124193119 Dzd = wuD = 1.00000000000001 0.9090909090

10、9091Vk= Vzd = wuV =0.99999999999993 0.90453403373329 0.904534033733350.99999999999995 0.30151134457776 0.301511344577781.00000000000000 -0.30151134457776 -0.30151134457776由输出结果可见，迭代次数k已经达到最大迭代次数max1=100，并且lambda的相邻两次迭代的误差Wc2.377 582,由wuV可以看出，lambda的特征向量Vk与真值Dzd的特征向量Vzd对应分量的比值相差较大，所以迭代序列发散.实际上，实数矩阵C的

11、特征值的近似值为，并且对应的特征向量的近似向量分别为=（0.90453403373329，0.30151134457776，-0.30151134457776），（-0.72547625011001,-0.21764287503300-0.07254762501100i, 0.58038100008801-0.29019050004400i），( -0.72547625011001, -0.21764287503300 + 0.07254762501100i, 0.58038100008801 + 0.29019050004400i), 是常数）.（4）输入MATLAB程序 D=-4 14 0

12、;-5 13 0;-1 0 2; V0=1,1,1; k,lambda,Vk,Wc=mifa(D,V0,0.00001,100), V,Dt =eig (D), Dtzd=max(diag(Dt), wuDt=abs(Dtzd-lambda), Vzd=V(:,2),wuV=V(:,2)./Vk,运行后屏幕显示结果请注意：迭代次数k,主特征值的近似值lambda,主特征向量的近似向量Vk,相邻两次迭代的误差Wc如下：k = lambda = Wc = 19 6.00000653949528 6.539523793591684e-006Dtzd = wuDt = 6.00000000000000

13、 6.539495284840768e-006Vk = Vzd = wuV =0.79740048053564 0.79740048053564 0.79740048053564 0.71428594783886 0.56957177181117 0.79740021980618 -0.24999918247180 -0.19935012013391 0.79740308813370（一） 原点位移反幂法的MATLAB主程序1用原点位移反幂法计算矩阵的特征值和对应的特征向量的MATLAB主程序1function k,lambdan,Vk,Wc=ydwyfmf(A,V0,jlamb,jd,max

14、1)n,n=size(A); A1=A-jlamb*eye(n); jd= jd*0.1;RA1=det(A1); if RA1=0disp(请注意：因为A-aE的n阶行列式hl等于零，所以A-aE不能进行LU分解.)returnendlambda=0;if RA1=0 for p=1:nh(p)=det(A1(1:p, 1:p);endhl=h(1:n);for i=1:nif h(1,i)=0disp(请注意：因为A-aE的r阶主子式等于零，所以A-aE不能进行LU分解.) returnendend if h(1,i)=0 disp(请注意：因为A-aE的各阶主子式都不等于零，所以A-aE

15、能进行LU分解.)k=1;Wc =1;state=1; Vk=V0;while(kjd)state=1;endk=k+1;%Vk=Vk2,mk=mk1,endif(Wc A=1 -1 0;-2 4 -2;0 -1 2;V0=1,1,1;k,lambda,Vk,Wc=ydwyfmf(A,V0,0.2,0.0001,10000)运行后屏幕显示结果请注意：因为A-aE的各阶主子式都不等于零，所以A-aE能进行LU分解.A-aE的秩R(A-aE)和各阶顺序主子式值hl、迭代次数k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相邻两次迭代的误差Wc如下：k = lambda = Wc =

16、 hl = 3 0.2384 1.0213e-007 0.8000 1.0400 0.2720Vk = V = D = 1.0000 -0.2424 -1.0000 -0.5707 5.1249 0 0 0.7616 1.0000 -0.7616 0.3633 0 0.2384 0 0.4323 -0.3200 -0.4323 1.0000 0 0 1.6367（2）输入MATLAB程序 A=1 -1;2 4;V0=20,1;k,lambda,Vk,Wc=ydwyfmf(A,V0,2.001,0.0001,100) 运行后屏幕显示结果请注意：因为A-aE的各阶主子式都不等于零，所以A-aE能进

17、行LU分解.A-aE的秩R(A-aE)和各阶顺序主子式值hl、迭代次数k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相邻两次迭代的误差Wc如下：k = lambda = Wc = hl = 2 2.0020 5.1528e-007 -1.0010 -0.0010Vk = V = D = 1.0000 -1.0000 0.5000 2 0 -1.0000 1.0000 -1.0000 0 3（3）输入MATLAB程序 A=-11 2 15;2 58 3;15 3 -3;V0=1,1,-1;k,lambdan,Vk,Wc=ydwyfmf(A,V0,8.26, 0.0001,100

18、) 运行后屏幕显示结果请注意：因为A-aE的各阶主子式都不等于零，所以A-aE能进行LU分解.A-aE的秩R(A-aE)和各阶顺序主子式值hl、迭代次数k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相邻两次迭代的误差Wc如下：k = lambdan= Wc = hl = 2 8.2640 6.9304e-008 -19.2600 -961.9924 -6.1256Vk = V = D = -0.7692 0.7928 0.6081 0.0416 -22.5249 0 0 0.0912 0.0030 -0.0721 0.9974 0 8.2640 0 -1.0000 -0.60

19、95 0.7906 0.0590 0 0 58.2609例 5.3.3 用原点位移反幂法的迭代公式（5.28），计算的分别对应于特征值，的特征向量，的近似向量，相邻迭代误差为0.001.将计算结果与精确特征向量比较.解 （1）计算特征值对应的特征向量的近似向量.输入MATLAB程序 A=0 11 -5;-2 17 -7;-4 26 -10;V0=1,1,1;k,lambda,Vk,Wc= ydwyfmf(A,V0,1.001, 0.001,100),V,D=eig(A);Dzd=min(diag(D), wuD= abs(Dzd- lambda),VD=V(:,1),wuV=V(:,1)./V

20、k,运行后屏幕显示结果请注意：因为A-aE的各阶主子式都不等于零，所以A-aE能进行LU分解.A-aE的秩R(A-aE)和各阶顺序主子式值hl、迭代次数k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相邻两次迭代的误差Wc如下：hl = -1.00100000000000 5.98500100000000 -0.00299600100000k = lambda = RA1 = 5 1.00200000000000 -0.00299600100000Vk = VD = wuV = -0.50000000000000 -0.40824829046386 0.816496580927

21、73 -0.50000000000000 -0.40824829046386 0.81649658092773 -1.00000000000000 -0.81649658092773 0.81649658092773Wc = Dzd = wuD = 1.378794763695562e-009 1.00000000000000 0.00200000000000 从输出的结果可见，迭代5次，特征向量的近似向量的相邻两次迭代的误差Wc1.379 e-009，由wuV可以看出，= Vk与VD 的对应分量的比值相等.特征值的近似值lambda 1.002与初始值1.001的绝对误差为0.001，而与的

22、绝对误差为0.002，其中,.（2）计算特征值对应特征向量的近似向量.输入MATLAB程序 A=0 11 -5;-2 17 -7;-4 26 -10;V0=1,1,1;k,lambda,Vk,Wc=ydwyfmf(A,V0,2.001, 0.001,100) ,V,D=eig(A); WD=lambda-D(2,2),VD=V(:,2),wuV=V(:,2)./Vk,运行后屏幕显示结果请注意：因为A-aE的各阶主子式都不等于零，所以A-aE能进行LU分解.A-aE的秩R(A-aE)和各阶顺序主子式值hl、迭代次数k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相邻两次迭代的误差

23、Wc如下：hl = -2.00100000000000 -8.01299900000000 0.00200099900000k = Wc = lambda = WD = 2 3.131363162302120e-007 2.00200000000016 0.00200000000016Vk = VD = wuV = -0.24999999999999 0.21821789023599 -0.87287156094401 -0.49999999999999 0.43643578047198 -0.87287156094398 -1.00000000000000 0.87287156094397

24、-0.87287156094397从输出的结果可见，迭代2次，特征向量的近似向量的相邻两次迭代的误差Wc3.131e-007，与的对应分量的比值近似相等.特征值的近似值lambda2.002与初始值2.001的绝对误差约为0.001，而lambda与的绝对误差约为0.002，其中,.（3）计算特征值对应特征向量的近似向量.输入MATLAB程序 A=0 11 -5;-2 17 -7;-4 26 -10;V0=1,1,1;k,lambda,Vk,Wc=ydwyfmf(A,V0,4.001, 0.001,100)V,D=eig(A); WD=lambda-max(diag(D),VD=V(:,3),

25、wuV=V(:,3)./Vk,运行后屏幕显示结果请注意：因为A-aE的各阶主子式都不等于零，所以A-aE能进行LU分解.A-aE的秩R(A-aE)和各阶顺序主子式值hl、迭代次数k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相邻两次迭代的误差Wc如下：hl =-4.00100000000000 -30.00899900000000 -0.00600500099999k = lambda = Wc = WD = 2 4.00199999999990 1.996084182914842e-007 0.00199999999990Vk = VD = wuV = 0.40000000

26、000001 -0.32444284226153 -0.81110710565380 0.60000000000001 -0.48666426339229 -0.81110710565381 1.00000000000000 -0.81110710565381 -0.81110710565381从输出的结果可见，迭代2次，特征向量的近似向量的相邻两次迭代的误差Wc1.996e-007，与的对应分量的比值近似相等.特征值的近似值的绝对误差近似为，而lambda与的绝对误差约为0.002，其中-0.400 000 000 000 00，-0.600 000 000 000 00，-1.000 00

27、0 000 000 00, .（二）原点位移反幂法的MATLAB主程序2用原点位移反幂法计算矩阵的特征值和对应的特征向量的MATLAB主程序2function k,lambdan,Vk,Wc=wfmifa1(A,V0,jlamb,jd,max1)n,n=size(A); jd= jd*0.1;A1=A-jlamb*eye(n);nA1=inv(A1); lambda1=0;k=1;Wc =1;state=1; U=V0;while(kjd)state=1;endk=k+1;endif(Wc A=0 11 -5;-2 17 -7;-4 26 -10;V0=1,1,1; k,lambda,Vk,W

28、c=wfmifa1(A,V0,1.001,0.001,100)运行后屏幕显示结果请注意迭代次数k,特征值的近似值lambda,对应的特征向量的近似向量Vk,相邻两次迭代的误差Wc如下：k = lambda = Wc =5 1.00200000000138 1.376344154436924e-006Vk = -0.50000000000000 -0.50000000000000 -1.00000000000000同理可得，另外与两个特征值对应的特征向量.（2）再用两种原点位移反幂法的MATLAB主程序，求对应的特征向量.输入MATLAB程序 A=0 11 -5;-2 17 -7;-4 26 -

29、10;V0=1,1,1;k,lambda,Vk,Wc=ydwyfmf(A,V0,0.99999999999997,0.001,100)运行后屏幕显示结果请注意：因为A-aE的各阶主子式都不等于零，所以A-aE能进行LU分解.A-aE的秩R(A-aE)和各阶顺序主子式值hl、迭代次数k,按模最小特征值的近似值lambda,特征向量的近似向量Vk,相邻两次迭代的误差Wc如下：hl = -0.99999999999997 6.00000000000045 0.00000000000010Vk = 0.50000000000000 0.50000000000000 1.00000000000000Wc

30、 = 4.317692037236759e-013RA1 = 1.039168751049192e-013k = 2lambda = 1.00000000000000输入MATLAB程序 A=0 11 -5;-2 17 -7;-4 26 -10;V0=1,1,1;k,lambda,Vk,Wc=wfmifa1(A,V0, 0.99999999999997,0.001,100)运行后屏幕显示结果Vk = 0.50000000000000 0.50000000000000 1.00000000000000请注意迭代次数k,特征值的近似值lambda,对应的特征向量的近似向量Vk,相邻两次迭代的误差W

31、c如下：k = 3lambda = 1.00000000000000Wc = 5.412337245047640e-0165.4 雅可比(Jacobi)方法及其MATLAB程序5.4.3 雅可比方法的MATLAB程序用雅可比方法计算对称矩阵的特征值和对应的特征向量的MATLAB主程序function k,Bk,V,D,Wc=jacobite(A,jd,max1)n,n=size(A);Vk=eye(n);Bk=A;state=1;k=0;P0=eye(n);Aij=abs(Bk-diag(diag(Bk);m1 i=max(Aij);m2 j=max(m1);i=i(j);while (kjd

32、)state=1;else returnendPk;Vk;Bk=B2;Wc;endif(kmax1)disp(请注意迭代次数k已经达到最大迭代次数max1,迭代次数k,对称矩阵Bk,以特征向量为列向量的矩阵V,特征值为对角元的对角矩阵D如下：)elsedisp(请注意迭代次数k,对称矩阵Bk,以特征向量为列向量的矩阵V,特征值为对角元的对角矩阵D如下：) endWc;k=k; V=Vk;Bk=B2;D=diag(diag(Bk);V1,D1 =eig(A,nobalance)5.常微分方程数值解法用matlab解微分方程组dx/dt=x-y-x(x2+y2)dy/dt=x+y-y(x2+y2)

33、x(0)=2y(0)1解析解：x,y=dsolve(Dx=x-y-x*(x2+y2), Dy=x+y-y*(x2+y2),x(0)=2,y(0)=1)得到的结果是解析解没有找到。用数值解。在Matlab下输入：edit，然后将下面两行百分号之间的内容，复制进去，保存 % function y=zhidao_rk4_5(t,x) %x,y变量分别用x(1),x(2)表示 y=x(1)-x(2)-x(1)*(x(1)2+x(2)2);x(1)+x(2)-x(2)*(x(1)2+x(2)2); % 在Matlab下面输入： t_end=10; x0=2;1; t,x=ode45(zhidao_rk4_5,0,t_end,x0); plot(t,x);legend(x,y);xlabel(t);figure;plot(x(:,1),x(:,2); xlabel(x); ylabel(y); 6.复化梯形公式用复化梯形公式求解 sinx积分，积分区间0pi建立Trapezoid.m文件function I,step= Trapezoid (f,a,b,eps) %f为函数，a为积分上限，b为积分下限，eps为积分精度，step为划分区间个数if(nargin=3) eps=1.0e-4;endn=1;h=(b-a)/2;I1=2;I