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

数值计算方法上机题目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)))<TOL
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);
p2 = r2+f*p1;
k = k + 1;
end
d = dot(r2,r2)/dot(p2,A*p2);
x = x+d*p2;
k = k + 1;
disp('方程组的conjgrad迭代次数为：')
k
disp('方程组的conjgrad迭代解为：')
x
disp('方程组的精确解为：')
xx=A\b
运行结果：
方程组的conjgrad迭代次数为：
k =
4
方程组的conjgrad迭代解为：
x =
1.00032203160110
1.00031075788355
1.00032203160110
方程组的精确解为：
xx =
1.00000000000000
1.00000000000000
1.00000000000000
比较：共轭梯度法迭代次数为4次，而G-S法的迭代次数需要54次，因此共轭梯度法计算量要小得多。