数值积分与线性方程组的解法

华北科技学院上机报告

系（部）建筑工程学院

专业、班级测绘B112 姓名学号 201105064226

课程名称数值分析

上机题目数值积分与线性方程组的解法

任课教师李慧

指导教师李慧

成绩（优、良、中、及格、不及格）

华北科技学院基础部

一．实验目的:

1）熟悉求解线性方程组以及数值积分的有关理论和方法；

2）会编制列主元消去法、LU分解法、平方根法、追赶法以及雅可比迭代和高斯-塞德尔迭代法的程序；

3）通过实际计算，进一步了解各种方法的优缺点，选择合适的数值方法，体会各种方法的精确度。

二．实验内容：

1．数值积分

梯形公式、辛普森公式、复化求积公式；

2．线性方程组求解

（1）高斯消去法、追赶法；

（2）雅可比迭代法、高斯塞德尔迭代法。

三、实验步骤与分析

1.数值积分的几种方法：

题目：已知积分精确值I=4.006994，分别用复化题型公式和复化辛普森公式计算其值。

I=⎰02)

+dx

1x

exp(

（1）.复化梯形公式

代码：

function I=trapez_v(f,h)

I=h*(sum(f)-(f(1)+f(length(f)))/2);

功能：复化求积公式进行函数积分

调用格式：I=trapez_v(f,h)

%f:等距节点上的函数值序列

%h:步长

程序如下：

clear

lexact=4.006994;

a=0;b=2;

fprintf('\n Extended Trapezoidal Rule\n');

fprintf(' n I Error\n');

n=1;

for k=1:6,n=2 * n;

h=(b-a)/n;i=1:n+1;

x=a+（i-1)*h; f=sqrt(1+exp(x));

I=trapez_v(f,h);

I=h*(sum(f)-(f(1）+f(length(f)))/2);

fprintf('%3.0f %10.5f %10.5f\n',n,I,Iexact-I);

end

结果:

Extended Trapezoidal Rule

n I Error

2 4.08358 -0.07659

4 4.02619 -0.01919

8 4.01180 -0.00480

16 4.00819 -0.00120

32 4.00729 -0.00030

64 4.00707 -0.00008

(2).复化辛普森公式

代码:

M文件：

function I=Simps_v(f,h)

n=length(f)-1;

if n==1,...

fprintf('Data has only one interval'),return; end

if n==2,...

I=h/3*(f(1)+4*f(2)+f(3));

return;end

I=0;

if n==3,...

I=3/8*h*(f(1)+3*f(2)+3*f(3)+f(4));

return;end

I=0;

if 2*floor(n/2)～=n,

I=3/8 * h * (f(n-2)+3*f(n-1)+3*f(n)+f(n+1));

m=n-3;

else

m=n;

end

I=I+(h/3)*(f(1)+4*sum(f(2:2:m))+f(m+1));

if m>2,I=I+(h/3)*2*sum(f(3:2:m));

end

function I=Simps_n(f_name,a,b,n)

h=(b-a)/n;

x=a+(0:n)*h;

f=feval(f_name,x);

I=Simps_v(f,h)

调用格式为: I=Simps_n('f_name',0,2,20)

结果为:

I=

4.0070

2．线性方程组的数值解法

（1）高斯消去法

题目：⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦

⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡-2313122434024121 ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡4321x x x x =⎥

⎥⎥⎥⎦

⎤⎢⎢⎢⎢⎣⎡6202813 代码:

M 文件：

function x = gauss(A, b)

n = length(b);

for k = 1 : n-1

if A(k,k)==0

fprintf('Error: the %dth pivot element equal to zero!\n',k); return;

end

index = [k+1:n];

m = -A(index,k)/A(k,k);

A(index,index) = A(index,index) + m*A(k,index);

b(index) = b(index) + m*b(k);

end

x = zeros(n,1);

x(n) = b(n)/A(n,n);

for i = n-1:-1:1

x(i) = ( b(i) - A(i,[i+1:n])*x([i+1:n]) )/A(i,i);

end

在Command Window 输入

>>A=[1 2 1 4;

2 0 4 3;

4 2 2 1;

-3 1 3 2];

b=[13,28,20,6]'

b =

13

28

20

6