数值积分与线性方程组的解法docx.docx
《数值积分与线性方程组的解法docx.docx》由会员分享，可在线阅读，更多相关《数值积分与线性方程组的解法docx.docx（11页珍藏版）》请在冰点文库上搜索。
数值积分与线性方程组的解法docx
华北科技学
上机报告
专业、班级测绘B112
姓名学号201105064226
课程名称数值分析
上机题目数值积分与线性方程组的解法
任课教师
指导教师李慧
成绩（优、良、中、及格、不及格）
华北科技学院基础部
1.实验目的:
1）熟悉求解线性方程组以及数值积分的有关理论和方法；
2）会编制列主元消去法、LU分解法、平方根法、追赶法以及雅可比迭代和高斯-塞徳尔迭代法的程序；
3）通过实际计算，进一步了解各种方法的优缺点，选择合适的数值方法，体会各种方法的精确度。
2.实验内容：
1.数值积分
梯形公式、辛普森公式、复化求积公式；
2.线性方程组求解
（1）高斯消去法、追赶法；
（2）雅可比迭代法、高斯塞徳尔迭代法。
三、实验步骤与分析
1.数值积分的几种方法：
题目：
已知积分精确值1=4.006994,分别用复化题型公式和复化辛普森公式计算其值。
[二J：
J]+exp（兀）dx
（1）•复化梯形公式
代码：
functionI二lrapez_v（f,h）
I=h*（sum（f）-（f（l）+f（length（f）））/2）;
功能：
复化求积公式进行函数积分
调用格式：
I=trapez_v（f,h）
%f:
等距节点上的函数值序列
%h:
步长
程序如下：
clear
lcxact=4.006994;
a=0;b=2;
fprintfC\nExtendedTrapezoidalRule\n‘）;
fprintf（'nIError\n，）;
n=l;
fork=l:
6,n二2*n;
h=（b-a）/n;i二1:
n+l;
x二a+（i-l）*h;f=sqrt（1+exp（x））;I=trapez_v（f,h）;
I=h*（sum（f）-（f
（1）+f（length（f）））/2）;fprintf（，%3.Of%10.5f%10.5f\n，,n,T,Texact-T）;end
结果：
Extended
TrapezoidalRule
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文件：
functionI=Simpsv（f,h）
n=length（f）T;
ifn二二1,...
fprintfCDatahasonlyoneintervaT）,return;
end
ifn==2,・・・
I=h/3*（f（l）+4*f
（2）+f（3））;
return;end
I二0;
ifn==3,…
I二3/8*h*（f（l）+3*f
（2）+3*f（3）+f（4））;
return;end
T=0;
if2*floor（n/2）〜二n,
1=3/8*h*（f（n-2）+3*f（n-l）+3*f（n）+f（n+l））;
m=n-3;
else
m=n;
end
1=1+（h/3）*（f
（1）+4*sum（f（2:
2:
m））+f（m+1））;
ifm>2,1=1+（h/3）*2*sum（f（3:
2:
m））;
endfunctionI二Simps_n（f_name,a,b,n）
h二（b-a）/n;
x二a+（0:
n）*h;
f=fcval（fname,x）;
I=Simps_v（f,h）
调用格式为:
I=Simps_n（，f_name，,0,2,20）结果为：
I二
4.0070
2.线性方程组的数值解法
（1）高斯消去法
1214
x\
13_
x2
28
x3
20
兀4
6
题目:
2043
4221
-3132
代码：
M文件:
functionx=gauss（A,b）n=length（b）;
fork=1:
n-1
ifA（k,k）==0
fprintf（JError:
the%dthpivotelementequaltozero!
\，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）;
fori=nT:
T:
1
x（i）=（b（i）—A（i,[i+l:
n]）*x（[i+l:
n]））/A（i,i）;
end
在CommandWindow输入
»A=[1214;
2043;
4221;
-3132];
b=[13,28,20,6]'
13
28
>>gauss（A,b）
结果：
ans=
3
-1
4
2
（2）追赶法
2-100
x\
~6_
x2
1
x3
-2
x4
_1_
-13-20
题目:
0-24-3
00-35
代码：
functionx=zhuiganfa
%首先说明：
追赶法是适用于三对角矩阵的线性方程组求解的方法，并不适用于其他类型矩阵。
%定义三对角矩阵A的各组成单元。
方程为Ax=d
%b为A的对角线元素仃~n）,a为-1对角线元素（2~n）,c为+1对角线元素
%A二[2-100
%-13-20
%0-24-3
%00-35]
a=[0-1-2-3];c=[-l-2-3];b二[2345];d=[61-21];n=length（b）;
u0=0;y0=0;a（l）=0;
%“追”的过程
L（l）=b（l）-a（l）*uO;
y（l）=（d（l）-yO*a（l））/L（l）;
u（l）二c（l）/L（l）;
fori=2:
（n-1）
L（i）=b（i）-a（i）*u（i-l）;
y（i）二（d（i）-y（i-l）*a（i））/L（i）;u（i）二c（i）/L（i）;
end
L（n）=b（n）-a（n）*u（n-1）;
y（n）=（d（n）-y（n-1）*a（n））/L（n）;
%“赶”的过程
x（n）=y（n）;
fori二（nT）：
-1:
1
x（i）二y（i）-u（i）*x（i+l）;
end
在命令窗口输入:
A=[2-1
0
0
■
-1
3
-2
0；
0
-2
4
-3;
0
0
-3
5];
a=[0-1-2-3];c二[-1-2-3];b二[2345];d二[61-21];
zhuiganfa
结果：
ans二
5.00004.00003.00002.0000
（3）雅克比迭代法和高斯赛德尔迭代法
题目：
分别用雅克比迭代法和高斯赛德尔迭代法求解线性方程组Ax二B,其中:
A二[4-11;4-8l;-215];B二[7-2115]'。
代码：
a.雅克比迭代
M文件：
functionX=jacobi（A,B,P,delta,maxi）
fprintf（，Tt.kx
（1）x
（2）x（3）err\n，）;N=length（B）;
fork=l:
maxi
forj=l:
N
x（j）=（B（j）-A（j,j+l:
N]）*P（[l:
j-l,j+l:
N]））/A（j,j）;
end
crr=abs（norm（x'-P））;
fprintf（'%3.Of,%10.6f,%10.6f,%10.6f,%10.6f\n，,k,x,err）;relerr=err/（norm（x）+eps）;
P二x';
if（errbreak
end
end
x二x'
命令窗口输入以下命令，并执行m文件。
»A二[4-11;4-81;-215];
B二[7-2115]';
P二[000]';
delta二0.0001;
maxl=12;
det（A）;
eig（A）;
X=jacobi（A,B,P,delta,maxi）;
计算结果：
Tt.
kx（l）
x
（2）x（3）
err
1,
1.750000,
2.625000,
3.000000,
4.353519
2,
1.656250,
3.875000,
3.175000,
1.265667
3,
1.925000,
3.850000,
2.887500,
0.394345
4,
1.990625,
3.948437,
3.000000,
0.163257
5,
1.987109,
3.995312,
3.006563,
0.047463
6,
1.997187,
3.994375,
2.995781,
0.014788
7,
1.999648,
3.998066,
3.000000,
0.006122
&
1.999517,
3.999824,
3.000246,
0.001780
9,
1.999895,
3.999789,
2.999842,
0.000555
10,
1.999987,
3.999927,
3.000000,
0.000230
2.0000
3.9999
3.0000
b.高斯赛德尔迭代
M文件：
functionX=jseid（A,B,P,delta,maxi）
N=length（B）
fork=l:
maxi
forj=l:
N
ifj==l
X（l）二（B
（1）-A（1,2:
N）*P（2:
N））/A（1,1）;
elseifj==N
X（N）=（B（N）-A（N,1:
N-1）*（X（1:
N-1））'）/A（N,N）;
else
X（j）二（B（j）-A（j,l:
j-l）*X（l:
j-l）-A（j,j+l:
N）*P（j+l:
N））/A（j,j）;
end
end
err=abs（norm（X'-P））;
relerr=err/（norm（X）+eps）;
P二X';
if（errbreak
end
end
X
命令窗口输入以下命令，并执行ni文件。
»A二[4-11;4-81;-215];
B=[7-2115]';
P=[000]';
delta二0.0000001;
max1=80;
X=gscid（A,B,P,delta,maxi）
计算结果:
X=
2.00004.00003.0000
四.实验总结
这次实验让我对MATLAB的使用更加熟练。
明白了如何定义M文件以及调用M文件得到结果。
也体会到数值微分的儿种方法各自的特点以及线性方程组儿种解法的优越性。
MATLAB软件在计算方而的确发挥了重要的作用，我应该对它有更进一步的学习，并熟练应用它來解决问题。
升级会员 