数值分析上机报告(2)

一．上机目的
1. 通过实际计算体会各种积分方法的精确度；会编写用龙贝格算法求定积分的程序。

2. 熟悉求解线性方程组的有关理论和方法；并会编制列主元消去法、LU 分解法。

二．上机环境
MATLAB 软件等。

三．上机内容
1.数值积分实验：复化积分、变步长积分；
2.线性代数方程组的直接解法：列主元消去法、LU 分解法。

四．实验内容
1、计算积分dx x ⎰+10211
①复化辛卜生公式
代码：
format long
fprintf('输入积分下限a 的值:\n')
a=input('a=');
fprintf('输入积分下上限b 的值:\n')
b=input('b=');
fprintf('输入区间分割次数n 的值:\n')
n=input('n=');
h=(b-a)/(2*n);
s1=0;s2=0;
i=a:h:b;
f=1./(1+i.^2);%被积函数
for t=1:n
s1=s1+f(2*t);
end
if n>=2
for t=1:n-1;
s2=s2+f(1+2*t);
end
else s2=0;
end
s=h/3*(f(1)+f(length(i))+4*s1+2*s2);
digits(6)
vpa(s)
结果：
输入积分下限a的值:
a=0
输入积分下上限b的值:
b=1
输入区间分割次数n的值:
n=2
ans =
0.785392
即：s(2)=0.785392
②变步长梯形法
代码：
a=0;b=1;%积分上下限
er=10^-6;%给定的误差
n=1;%初始分割数
h=b-a;%初始步长
T2=h/2*(1/(1+a^2)+1/(1+b^2)); T1=0;%给T1赋初值，保证进入循环while abs(T2-T1)>=er
s=0;
T1=T2;
for k=0:n-1
x=a+k*h+h/2;
s=s+1/(1+x^2);
T2=1/2*(T1+h*s);
end
h=h/2;
n=2*n;
end
digits(6)
fprintf('输出分割段数：')
disp(n)
fprintf('输出积分值：')
j=vpa(T2)
结果：
输出分割段数： 512
输出积分值：
j =
0.785398
2、 用列主元消去法求解方程组Ax=b
⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡-----=38126
5973274581221A ,b=⎥⎥⎥⎥⎦
⎤⎢⎢⎢⎢⎣⎡4911427。

代码： format short
%列主元高斯消去法
A=[1 2 -12 8;5 4 7 -2;-3 7 9 5;6 -12 -8 3];%系数矩阵 B=[27;4;11;49];%方程右端值
%选主元
n=length(A);%矩阵维数
for r=1:n
a=0;
for i=r:n
if a<=abs(A(i,r))
a=abs(A(i,r));
t=i;
end
end
%换行
if t~=r
p(1,:)=A(r,:);
A(r,:)=A(t,:);
A(t,:)=p(1,:);
q=B(r);
B(r)=B(t);
B(t)=q;
end
%消元
for i=r:n-1
for j=i+1:n
m(j)=A(j,r)/A(i,i);
A(j,:)=-A(i,:)*m(j)+A(j,:); B(j)=-m(j)*B(i)+B(j);
end
end
end
%求解x
for i=1:n
t=n-i+1;
x(t)=B(t);
for t1=t+1:n
x(t)=x(t)-A(t,t1)*x(t1);
end
x(t)=x(t)/A(t,t);
end
fprintf('输出x 解矩阵：\n')
disp(x')
结果：
输出x 解矩阵：
3.0000
-2.0000
1.0000
5.0000
3、用LU 分解法求解方程组Ax=b
⎥⎥⎥⎦
⎤⎢⎢⎢⎣⎡--=122140
111
A ⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=156b 代码：
A=[1 1 1;0 4 -1;2 -2 1];%系数矩阵 B=[6;5;1];%方程右端值
n=length(A);%矩阵维数
U(1,:)=A(1,:);%U 矩阵第一行
for i=2:n
L(i,1)=A(i,1)/A(1,1);%L 矩阵第一列 end
%求出L,U 矩阵
for i=1:n
L(i,i)=1;
end
sum=0;sum1=0;
for r=2:n
sum=0;
for i=r:n
for k=1:r-1
sum=sum+L(r,k)*U(k,i);
end
U(r,i)=A(r,i)-sum;
end
if r~=n
sum1=0;
for i=r+1:n
for k=1:r-1
sum1=sum1+L(i,k)*U(k,r);
end
L(i,r)=(A(i,r)-sum1)/U(r,r);
end
end
end
%计算y的值
y(1)=B(1);
for i=2:n
sum2=0;
for k=1:i-1
sum2=sum2+L(i,k)*y(k);
end
y(i)=B(i)-sum2;
end
%计算x的值
x(n)=y(n)/U(n,n);
for t=1:n-1
i=n-t;
sum3=0;
for k=i+1:n
sum3=sum3+U(i,k)*x(k);
end
x(i)=(y(i)-sum3)/U(i,i);
end
fprintf('输出L矩阵：\n')
disp(L)
fprintf('输出U矩阵：\n')
disp(U)
fprintf('输出x解矩阵：\n')
disp(x')
结果：
输出L矩阵：
1 0 0
0 1 0
2 -1 1
输出U矩阵：
1 1 1
0 4 -1
0 0 -2
输出x解矩阵：
1
2
3。