计算方法实验.docx

1、计算方法实验计算方法实验报告 目录计算方法实验报告 - 1 -实验一 舍入误差与数值稳定性 - 2 -实验二 方程求根 - 3 -实验三 线性方程组数值解法 - 4 -实验四 插值法 - 8 -实验五 曲线拟合 - 9 -实验六 数值积分与数值微分 - 11 -实验一 舍入误差与数值稳定性设=已知其精确值为(.1）编制从大到小的顺序计算的程序；2）编制按从小到大的顺序计算的程序；3）按2种顺序分别计算，并指出有效位数。思路分析：两种不同的累加方式。解：1)程序如下：（预设N=100）N=100;S=0;while(N-1) S=S+1/(N2-1); N=N-1;endfprintf(%f,S

2、);2）程序如下： N=100;S=0;for j=2:N S=S+1/(j2-1);endfprintf(%f,S);3）降序：= 0.749000 = 0.749900 =0.749967升序：=0.749000 =0.749900 =0.749967实验二 方程求根用牛顿法求下列方程的根：思路分析：首先应是牛顿迭代式，之后判断收敛即可。解：程序如下：eps=5e-6;delta=1e-6;N=100;k=0;x0=1;while(1) x1=x0-func1(x0)/func2(x0); k=k+1; if(kN|abs(x1)eps) disp(newton method failed

3、) break; end if abs(x1)1 d=x1-x0; else d=(x1-x0)/x1; end x0=x1; if(abs(d)eps|abs(func1(x1)delta) break; endendfprintf(%f,x0);function y=func2(x) y=2*x-exp(x);endfunction y=func1(x) y=x*x-exp(x);end运行结果：-0.703467实验三 线性方程组数值解法1. 编写用追赶法解三对角线性方程组的程序，并解下列方程：Ax=b,其中，=-4 1 =-271 -4 1 -15 1 -4 1 -15 1 -4 1

4、-15 1 -4 -15思路分析：按照追赶法的递推公式求解。 解：程序如下: function x=Trid(a,b,c,d)% 追赶法求解三对角的线性方程组 Ax=d% b为主对角线元素，a，c分别为次对角线元素，d为右端项 % A= b1 c1% a2 b2 c2% .% a_(n-1) b_(n-1) c_(n-1)% a_(n) b_(n) % b=b1.b_(n)% a=0 a2.a_(n)% c=c1.c_(n-1)n=length(b);u(1)=b(1);for i=2:n l(i)=a(i)/u(i-1); u(i)=b(i)-l(i)*c(i-1);endy(1)=d(1)

5、;for i=2:n y(i)=d(i)-l(i)*y(i-1);endx(n)=y(n)/u(n);for i=n-1:-1:1 x(i)=(y(i)-c(i)*x(i+1)/u(i);endfor i=1:n fprintf(x%d=%fn,i,x(i);end运行结果：x1=8.705758, x2=7.823032, x3=7.586371, x4=7.522453, x5=7.503440,x6=7.491306,x7=7.461785,x8=7.355835,x9=6.961556, x10=5.4903892. 分别用雅克比迭代法与高斯-赛德尔迭代法解下列方程组：RI=V,其中R

6、=31 -13 0 0 0 -10 0 0 0 V=-15 -13 35 -9 0 -11 0 0 0 0 27 0 -9 31 -10 0 0 0 0 0 -23 0 0 -10 79 -30 0 0 0 -9 0 0 0 0 -30 57 -7 0 -5 0 -200 0 0 0 7 47 -30 0 0 12 0 0 0 0 0 -30 41 0 0 -7 0 0 0 0 -5 0 0 27 -2 7 0 0 0 0 0 0 0 -2 29 -10思路分析：分别按照雅克比与高斯迭代法递推式编制程序。解：程序如下：雅克比迭代法function jacobi(A,b,max,eps)n=le

7、ngth(A);x=zeros(n,1);x1=zeros(n,1);k=0;while 1 x1(1)=(b(1)-A(1,2:n)*x(2:n,1)/A(1,1); for i=2:n-1 x1(i)=(b(i)-A(i,1:i-1)*x(1:i-1,1)-A(i,i+1:n)*x(i+1:n,1)/A(i,i); end x1(n)=(b(n)-A(n,1:n-1)*x(1:n-1,1)/A(n,n); k=k+1; if sum(abs(x1-x)=max fprintf(the method is disconvergentn); break; end x=x1;endif kmax

8、 for i=1:n fprintf(x%d=%fn,i,x1(i); endend 运行结果： x1=-0.200550,x2=0.368393,x3=-0.731860,x4=-0.300318, x5=-0.446577,x6=0.399384,x7=0.121500,x8=0.151792,x9=-0.334359高斯-赛德尔迭代法function GauseSeidel(A,b,max,eps)n=length(A);x=zeros(n,1);x1=zeros(n,1);k=0;while 1 x1(1)=(b(1)-A(1,2:n)*x(2:n,1)/A(1,1); for i=2

9、:n-1 x1(i)=(b(i)-A(i,1:i-1)*x1(1:i-1,1)-A(i,i+1:n)*x(i+1:n,1)/A(I,I); end x1(n)=(b(n)-A(n,1:n-1)*x1(1:n-1,1)/A(n,n); k=k+1; if sum(abs(x1-x)=max fprintf(the method is discomvergentn); break; end x=x1;endif kmax for i=1:n fprint(x%d=%fn,i,x1(i); endend运行结果：x1=-0.200550,x2=0.368393,x3=-0.731860,x4=-0.

10、300318, x5=-0.446577,x6=0.399384,x7=0.121500,x8=0.151792,x9=-0.334359实验四 插值法按下列数据做五次插值，并求x1=0.46,x2=0.55,x3=0.60时的函数近似值。Xi0.300.420.500.580.660.72yi1.044031.084621.118031.156031.198171.23223思路分析：用牛顿插值法多项式来解决此问题。解：程序如下：function newtint(x,y,xhat)n=length(y);c=y(:);for j=2:n for i=n:-1:j c(i)=(c(i)-c(i

11、-1)/(x(i)-x(i-j+1); endendyhat=c(n);for i=n-1:-1:1 yhat=yhat*(xhat-x(i)+c(i);endfprintf(N(%f)=%f,xhat,yhat);运行结果：x=0.46，y=1.100742;x=0.55,y=1.141271;x=0.60,y=1.166194;实验五 曲线拟合试分别用抛物线y=a+bx+c和指数曲线y=a拟合下列数据：X11.522.533.544.5y33.479.5122.65159.05189.15214.15238.65252.50x55.566.577.58y267.55280.5296.653

12、01.40310.4318.15325.15比较2个拟合函数的优劣。解： 1.抛物线拟合 程序如下： function ZXE(x,y,m)S=zeros(1,2*m+1);T=zeros(m+1,1);for k=1:2*m+1 S(k)=sum(x.(k-1);endfor k=1:m+1 T(k)=sum(x.(k-1).*y);endA=zeros(m+1,m+1);a=zeros(m+1,1);for i=1:m+1 for j=1:m+1 A(i,j)=S(i+j-1); endenda=AT;for k=1:m+1 fprintf(a%d=%fn,k,a(k); end运行结果：

13、a1=-45.333297，a2=94.230200，a3=-6.131610；即：y=-45.333297+94.230200x-6.1316102.指数曲线拟合 程序如下：function ZXEI(x,y,m)y=log(y);S=zeros(1,2*m-1);T=zeros(m,1);for k=1:2*m-1 S(k)=sum(x.(k-1);endfor k=1:m T(k)=sum(x.(k-1).*y);endA=zeros(m,m);a=zeros(m,1);for i=1:m for j=1:m A(i,j)=S(i+j-1); endenda=AT;a(1)=exp(a(

14、1);for k=1:m fprintf(a%d=%fn,k,a(k);end运行结果：a1=67.402,a2=0.238960，即：y=67.402实验六 数值积分与数值微分用复化梯形公式和复化辛普森公式计算积分和观察n为多少时所得近似值具有6位有效数字。解： MATLAB程序如下： 复化梯形公式function trap(a,b)n=1030;for i=1:3 x=a:(b-a)/n:b; h=(b-a)/n; s=(h/2)*(f(a)+2*sum(f(x(2:1:n-1)+f(b); fprintf(s(%d)=%fn,n,s); n=n+1;endfunction y=f(x)y

15、=(1+cos(x).2).0.5;运行结果：s(1030)=1.908574s(1031)=1.908575s(1032)=1.908577n=1030时才有6位有效数字。 复化辛普森公式function simpson(a,b)n=2;for i=1:3 x=a:(b-a)/(2*n):b; m=2*n+1; h=(b-a)/n; s=(h/6)*(f(a)+2*sum(f(x(3:2:m-2)+4*sum(f(x(2:2:m-1)+f(b); fprintf(s(%d)=%fn,n,s); n=n*2;endfunction y=f(x)y=(1+cos(x).2).0.5;运行结果：s

16、(2)=1.910141s(4)=1.910099s(8)=1.910099 n=4时即有6位有效数字； 复化梯形公式function trap(a,b)n=1954;for i=1:3 x=a:(b-a)/n:b; h=(b-a)/n; s=(h/2)*(1+2*sum(f(x(2:1:n-1)+f(b); fprintf(s(%d)=%fn,n,s); n=n+1;endfunction y=f(x)y=tan(x)/x;运行结果：s(1954)=0.848456s(1955)=0.848456s(1956)=0.848456n=1954时才有6位有效数字。 复化辛普森公式function simpson(a,b)n=4;for i=1:3 x=a:(b-a)/(2*n):b; m=2*n+1; h=(b-a)/n; s=(h/6)*(1+2*sum(f(x(3:2:m-2)+4*sum(f(x(2:2:m-1)+f(b); fprintf(s(%d)=%fn,n,s); n=n*2;endfunction y=f(x)y=tan(x)./x;运行结果：s(4)=0.848972s(8)=0.848967s(16)=0.848967 n=8时即有6位有效数字；