数值积分matlab实验程序.docx

1、数值积分matlab实验程序这是三次样条插值，及龙格现象的图。这里版权申明：程序版权归崛起强（这是XXID）所有。其中只有部分程序摘自精通matlab科学计算，但均作出了注解及修改。请勿转载！请勿用于商业用途！用所编写的程序计算积分， 复合梯形公式、复合辛普森公式、龙贝格公式、高斯列让德公式、高斯拉盖尔求积公式function I,n=RecursionTrapezoid(f,a,b,error)%n为节点数if nargin3 disp(You must input three variables at least); %quit;用quit不好，直接退出matlab，改为return,才能

2、看到disp。 return;elseif nargin=3 error=1.0e-4;%error是个精度eps,精度太高难以等待1.0e-10相当久endh=b-a;%步长%y0=0.5*h*(fun(a)+fun(b);%一步%n=1; %增加的节点数%flag=1;%while(flag)% sum=0;% xs=a+0.5*h;% for i=1:n;%默认的做了步插值点。翻倍插值% sum=sum+fun(xs);% xs=xs+h;% end % y1=0.5*y0+0.5*h*sum;%y0是对前面结果的复用，sum是新节点之和，新步长0.5h% if abs(y1-y0)er

3、ror sum=0; xs=a+0.5*h; y0=y1;%后推 for i=1:n sum=sum+f(xs); xs=xs+h; end y1=0.5*y1+0.5*h*sum;%1式 n=2*n;%保证循环变量才置后，且不引起不一致性 h=h/2; %y1=y0; %如果与1式合并，无法控制循环endI=y1;End计算上三式得值：I,n=RecursionTrapezoid(x0.5*log(x),0,1)I = NaNn = 1都无法求解，前两个是log(0) ，sin(0)/0无法计算，后者是无穷积分。采用折中的方法是，二者都在积分区间一致有界，误差可估计，故用1.0e-6代替0求

4、近似解运行都超过十分钟，终止计算。采用内联函数inline方法加快效率（代替sym）第二式I = 0.946057561038914 n = 32 I = 0.946082070343826 n = 32768（为1.0e-10）I = -0.444395950087616 n = 1024 第一式n为插值点数复合辛普森公式function I,step=IntSimpson(f,a,b,type,eps)%辛普森系列公式求函数f在区间a,b上的定积分if(type=3 &nargin=4) disp(lack varargin);endI=0;switch type case 1,%辛普森公

5、式 I=(b-a)/6)*(f(a)+. 4*f(a+b)/2)+. f(b); step=1; case 2,%辛普森3/8公式 I=(b-a)/8)*(f(a)+. 3*f(2*a+b)/3)+. 3*f(a+2*b)/3)+f(b); step=1; case 3,%复合辛普森公式 n=2; h=(b-a)/2; I1=0; I2=(f(a)+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)*(f(x)+. 4*f(x+x1)/2)+. f

6、(x1); endendI=I2;step=n;endI,step=IntSimpson(f,a,b,type,eps)计算前二式得：I = -0.444332482271752 n = 164（n为2分次数） I = 0.946082310887237 n = 4 龙贝格公式function I,step=Romberg(f,a,b,error)if nargin3 disp(You must input three variables at least); return;elseif nargin=3 error=1.0e-4;%1.0e-10不现实end%h=b-a;%preT(1)=0

7、.5*h*(f(a)+f(b);%可扩增向量%exitflag=0;%k=2;%n=1;%增加的节点数%while(exitflag)%该算法则指用到上两行：对角行和次对角行 % xs=a+0.5*h;% sum=0.;% for i=1:n% sum=sum+f(xs);% xs=xs+h;% end% postT(1)=0.5*preT(1)+0.5*h*sum;% for i=2:k% postT(i)=4(i-1)/(4(i-1)-1)*postT(i-1)-1/(4(i-1)-1)*preT(i-1);% end% if(abs(postT(k)-preT(k-1)error k=k

8、+1; h=h/2; Q=0;%sum for i=1:M x=a+h*(2*i-1); Q=Q+f(x); end T(k+1,1)=T(k,1)/2+h*Q;%梯形步长减半的复合。这是第二步 M=2*M;%节点数翻倍 for j=1:k %这是把整个下三角都求出来。 T(k+1,j+1)=T(k+1,j)+(T(k+1,j)-T(k,j)/(4j-1);%简单的变换，分离个1出 %来，输入更方便，下标从1开始故j+1 end tol=abs(T(k+1,j+1)-T(k,j);endI=T(k+1,k+1);step=k;%第一列k次二分。EndI = -0.444444214475247

9、 step = 14 I = 0.946082070387222 step = 3高斯拉盖尔求积公式function I=IntGaussLager(f,n,AK,XK)%高斯拉盖尔，n大于5自加精度,AK,XK,可自定义，默认最多是五项。if(n6 & nargin=2) AK=0; XK=0;else I=sum(AK.*f(XK);endswitch n case 2, I=0.853553*f(-0.585786)+. 0.146447*f(3.414214); case 3, I=0.711093*f(0.415575)+. 0.278518*f(2.294280)+. 0.0103

10、893*f(6.289945); case 4 I=0.603154*f(0.322548)+. 0.357419*f(1.745761)+. 0.0388879*f(4.536620)+. 0.000539295*f(9.395071); case 5 I=0.521756*f(0.263560)+. 0.398667*f(1.413403)+. 0.0759424*f(3.596426)+. 0.00361176*f(7.085810)+. 0.0000233700*f(12.640801);End求解第三式：ans = 0.498903451386692高斯列让德公式function I

11、=IntGauss(f,a,b,n,AK,XK)%AKXK1if(n load Ag.dat AgAg = 10.000000000000000 -7.000000000000000 0 1.000000000000000 -3.000000000000000 2.099999000000000 6.000000000000000 2.000000000000000 5.000000000000000 -1.000000000000000 5.000000000000000 -1.000000000000000 2.000000000000000 1.000000000000000 0 2.0

12、00000000000000 b1=8 5.900001 5 1b1 = 8.000000000000000 5.900001000000000 5.000000000000000 1.000000000000000 GaussOrder(Ag,b1)ans = 1.0e+007 * 0.000000800000000 0.000000830000100 2.075000349709961 0.000000100000000 Gauss_pivot(Ag,b1)ans = 0.000000000000000 -1.000000000000000 1.000000000000000 1.0000

13、00000000000差异：显然用高斯顺序消去法得到的结果是错误的，用高斯列主元消去法得到了正确解。造成这种区别的原因可能是：后者避免了微小量作为列主元（即分母），不会放大误差，保证了算法的稳定性，避免了误差危害高斯顺序消元法function x=GaussOrder(A,b)if nargin=2 A=A,b;endN=size(A);s=N(1);if N(1)=N(2)-1&length(N)=2 disp(error-input); return;endfor i=1:s %A约化为上三角形矩阵 if A(i,i)=0 disp(对角元素有0); return; end for j=i

14、+1:s A(j,i)=-A(j,i)/A(i,i); for k=i+1:N(2) A(j,k)=A(j,k)+A(j,i)*A(i,k); end endend function X=BackSubstitution(A) %这里用U，b不方便 N=size(A); s=N(1); X=ones(s,1); A(s,s+1)=A(s,s+1)/A(s,s); for i=s-1:1 for j=s:i+1 A(i,s+1)=A(i,s+1)-A(j,s+1)*A(i,j); end A(i,s+1)=A(i,s+1)/A(i,i); end for i=1:s X(i)=A(i,s+1);

15、 end endx=BackSubstitution(A);end列主元消去法function X,det=Gauss_pivot(A,b)if nargin=0 disp(没参数) return;endN=size(A);n=length(b);if N(1)=N(2)&length(N)=2 disp(输入A有误); return;endif N(1)=n disp(A和b不匹配); return;enddet=1;max=0;t=0;c=0;c2=0;for i=1:n %选列主元 for j=i:n%找i列主元 if maxabs(A(j,i) max=A(j,i); t=j; end

16、 end if max load Ao1.dat b2=1 1 1b2 = 1 1 1 Gauss_pivot(Ao1,b2)ans = 1.0e+003 * 1.592599624841381 -0.631911376202549 -0.493617724759390 load Ao2.dat Gauss_pivot(Ao2,b2)ans = 1.0e+002 * 1.195273381259593 -0.471426044312964 -0.368402561091259差异：相差一个数量级，可见微小的系数差别可能造成结果的巨大差异，这是因为消元过程中不同量得权值在变化，且变化巨大，每次消

17、元都有误差积累。3. 编写追赶法的程序代码，并求解教材第177页的第9题。 load flp.dat followup(flg,b3)ans = 0.833333333333333 0.666666666666667 0.500000000000000 0.333333333333333 0.166*6667雅可比迭代法function y=Jacobi(A,b)m,n=size(A);if m=n disp(error); quit;endif m=length(b) disp(error); quit;endx0=zeros(n,1);flag=1;while flag for i=1:n

18、 x1(i)=(b(i)-A(i,:)*x0)/A(i,i)+x0(i); end if norm(x1.-x0.)1.0e-7 flag=0; end x1=x0;end y=x1;end高斯塞德尔function GaussSeidel(A,b,x0,eps,M)if nargin=3 eps=1.0e-6; M=200;elseif nargin=4 M=200;elseif nargin=eps x0=x; x=G*x0+f; n=n+1; if(n=M) disp(Warning:迭代次数太多，可能不收敛); return; endendSOR迭代法function x,n=SOR(

19、A,b,x0,w,eps,M)if nargin=4 eps=1.0e-6; M=200;elseif nargin4 disp(error); return;elseif nargin=5 M=200;endif(w=2) disp(error); return;endD=diag(diag(A);L=-tril(A,-1);U=-triu(A,1);B=inv(D-L*w)*(1-w)*D+w*U);f=w*inv(D-L*w)*b;x=B*x0+f;n=1;while norm(x-x0)=eps x0=x; x=B*x0+f; n=n+1; if(n=M) disp(warning:迭

20、代次数太多，可能不收敛); return; endend1. 编写雅克比迭代法、高斯-赛德尔迭代法的程序代码；取相同的初始点，分别用所得的程序求解如下方程组 比较两种方法所得的结果；如果有差异，解释结果差异的原因。 ans =（初值（0,-1,2） 雅可比迭代的两个结果 -0.135*5135 ans = -3 3 1（初值（0 0 0）.0810* 3.918918918918919高斯-塞德尔：1. x = 2. x = -0.135*5135 -3.0810* 3 3.918918918918919 1 2. 用SOR法的程序解如下线性方程组，松弛因子分别为要求当时迭代终止。比较不同采用松弛因子时的迭代次数。ans =（w取0.9且初值是0向量） ans =（w取1.0且初值是0向量） -4.000000039257963 - 4.000000231430727 3.000000001923201 2.999999856370402 2.00