数值分析作业复化求积公式

数值计算方法上机题目3

计算定积分的近似值：

2

2

1x e xe dx =⎰ 要求：

（1）若用复化梯形公式和复化Simpson 公式计算，要求误差限7102

1-⨯=ε，分别利用他们的余项估计对每种算法做出步长的事前估计；

（2）分别利用复化梯形公式和复化Simpson 公式计算定积分；

（3）将计算结果与精确解比较，并比较两种算法的计算量。

解：

（1) x xe x f =)(,所以x x k xe ke x f +=)()(，x x xe e x f +=2)(''，x x xe e x f +=4)()4( x x xe e x f +=6)()6(

对于复化梯形公式： )(12)(''2ηf h a b f R n --=，2max ''4)(e f =η，n

h 1= 代入数据可知 722102

1124-⨯≤n e ，57.7018≥n 取7019=n

对于复化Simpson 公式 )()2(180)()4(4ηf h a b f R n --=，2max )4(6)(e f =η，n

h 1= 代入数据可知 742102

128806-⨯≤n e ，56.23≥n 取24=n

（2）复化梯形公式：

函数

function y=fun(x)

y=x*exp(x);

程序：

clc

Clear

% 复化梯形计算

format long

a=1;b=2;

n=7019;% 区间划分为m等份

h=(b-a)/n;% 步长，根据误差限由该算法的余项作事前估计得到

ty1=fun(a)+fun(b);

ty2=0;

for i=1:n-1

x=a+i*h1;

ty2=ty2+fun(x);

end

T=h*(ty1+2*ty2)/2;

T

对于复化Simpson公式

clc

Clear

% 复化Simpon计算

format long

a=1;b=2;

n=24;% 区间划分为n等份

h=(b-a)/(2*n);% 步长，根据误差限由该算法的余项作事前估计得到sy1=fun(a)+fun(b);

sy2=0;sy3=0;

for j=1:2*n-1

x=a+j*h2;

if rem(j,2)==0

sy3=sy3+fun(x);

else

sy2=sy2+fun(x);

end

end

S=h*(sy1+4*sy2+2*sy3)/3;

S

% 精确值

Exactanswer=exp(2)

运算结果：

T =

7.389056127230221

S =

7.389056126214707

Exactanswer =

7.389056098930650

（3）比较可知复化梯形公式的计算量较大

1．用共轭梯度法和G-S 迭代法分别求解下面的方程组：

⎥⎥⎥⎦

⎤⎢⎢⎢⎣⎡-=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-14514103131021310321x x x 迭代20次或满足()(1)11

10k k x x --∞-<时停止计算。

解：

（1）G-S 迭代法程序：

% G-S 迭代计算

clc

clear

A=[10 3 1;2 -10 3;1 3 10];

b=[14 -5 14]';

x0=[0 0 0]';

TOL=0.00000000001;

n=length(b);

k=0;

x=zeros(n,1);

while 1

x(1)=(b(1)-A(1,2:n)*x0(2:n))/A(1,1);

for i=2:n-1

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

end

x(n)=(b(n)-A(n,1:n-1)*x(1:n-1))/A(n,n);

k=k+1;

if max(abs(x(1:n)-x0(1:n)))

break;

end

x0=x;

end

disp('方程组的G-S 迭代次数为：')

k

disp('方程组的G-S 迭代解为：')

x

disp('方程组的精确解为：')

xx=A\b

运行结果：

方程组的G-S 迭代次数为：

k =

17

方程组的G-S迭代解为：

x =

0.99999999999970

0.99999999999978

1.00000000000009

方程组的精确解为：

xx =

1.00000000000000

1.00000000000000

1.00000000000000

（2）共轭梯度法程序：

% 共轭梯度法计算

clc

clear

A=[10 3 1;2 -10 3;1 3 10];

b=[14 -5 14]';

x0=[0 0 0]';

if(nargin == 3)

eps =1.0e-11;

end

r1=b-A*x0;

p1=r1;

d=dot(r1,r1)/dot(p1,A*p1);

x=x0+d*p1;

r2=r1-d*A*p1;

f=dot(r2,r2)/dot(r1,r1);

p2=r2+f*p1;

k=1;

for(i=1:(rank(A)-1))

x0 = x;

p1 = p2;

r1 = r2;

d = dot(r1,r1)/dot(p1,A*p1);

x = x0+d*p1;

r2 = r1-d*A*p1;

f = dot(r2,r2)/dot(r1,r1);