物探数据反演上机实验.docx
《物探数据反演上机实验.docx》由会员分享，可在线阅读，更多相关《物探数据反演上机实验.docx（27页珍藏版）》请在冰豆网上搜索。
物探数据反演上机实验
中南大学实验报告
物探数据处理与反演上机实验
班级：
姓名：
学号：
指导老师：
2015.6.27
物探数据处理上机实验
一．插值
1.拉格朗日插值
实例：
根据下面的数据点求出其拉格朗日插值格式，并计算当x=1.6时y的值。
x
0
0.5
1.0
1.5
2.0
2.5
3.0
y
0
0.4794
0.8145
0.9975
0.9093
0.5985
0.1411
function[f,f0]=Language（x,y,x0）
%求已知数据点的拉格朗日插值多项式
%已知数据点的x坐标向量：
x
%已知数据点的y坐标向量：
y
%插值点的x坐标：
x0
%求得的拉格朗日插值多项式：
f
%x0处的插值：
f0
symst;
if（length（x）==length（y））
n=length（x）;
else
disp（'x和y的维数不相等！
'）;
return;
end%检错
f=0.0;
for（i=1:
n）
l=y（i）;
for（j=1:
i-1）
l=l*（t-x（j））/（x（i）-x（j））;
end;
for（j=i+1:
n）
l=l*（t-x（j））/（x（i）-x（j））;%计算拉格朗日基函数
end;
f=f+l;%计算拉格朗日插值函数
simplify（f）;%化简
end
f0=subs（f,'t',x0）;%计算插值点的函数值
运行程序；
x=0:
0.5:
3;
y=[00.47940.84150.99750.90930.59850.1411];
[f,f0]=Language（x,y,1.6）%计算输出的拉格朗日插值多项式
%计算x=1.6时的插值输出值f0
运行结果
f=
-799/3125*t*（t-1）*（t-3/2）*（t-2）*（t-5/2）*（t-3）+561/500*t*（t-1/2）*（t-3/2）*（t-2）*（t-5/2）*（t-3）-133/75*t*（t-1/2）*（t-1）*（t-2）*（t-5/2）*（t-3）+3031/2500*t*（t-1/2）*（t-1）*（t-3/2）*（t-5/2）*（t-3）-399/1250*t*（t-1/2）*（t-1）*（t-3/2）*（t-2）*（t-3）+1411/112500*t*（t-1/2）*（t-1）*（t-3/2）*（t-2）*（t-5/2）
f0=
0.9996
2.两次样条插值
实例：
根据下面的数据点求出二次样条插值多项式，并计算x=4.3时y的值。
x
1
2
3
4
5
6
7
8
y
0.8415
0.9093
0.1411
-0.7568
-0.9589
-0.2794
0.6570
0.9894
y’
0.5403
-0.4161
-0.9900
-0.6536
0.2837
0.9602
0.7539
-0.1455
function[f,f0,fd0]=SecSample（x,y,y_1,x0）
symst;
f=0.0;
f0=0.0;
if（length（x）==length（y））
if（length（y）==length（y_1））
n=length（x）;
else
disp（'y和y的倒数的维数不相等！
'）;
return;
end
else
disp（'x和y的倒数的维数不相等！
'）;
return;
end
fori=1:
n
if（x（i）<=x0）&&（x（i+1）>=x0）
index=i;
break;
end
end
d=y_1
（1）*（x
（2）-x
（1））/2+y
（1）;
fori=2:
n-1
d=d+y_1（i）*（x（i+1）-x（i-1））/2;
end
h=x（index+1）-x（index）;
f=y_1（index+1）*（t-x（index））^2/2/h+...
y_1（index）*（t-x（index+1））^2/2/h+d;
fd=（t-x（index））*y_1（index+1）/h+y_1（index）*（x（index+1）-t）/h;
f0=subs（f,'t',x0）;
fd0=subs（fd,'t',x0）;
end
x=1:
8;
y=[0.84150.90930.1411-0.7568-0.9589-0.27940.65700.9894]
y_1=[0.5403-0.4161-0.9900-0.65360.28370.96020.7539-0.1455]
[f,f0,fd0]=SecSample（x,y,y_1,4.3）
运行结果
f=
2837/20000*（t-4）^2-817/2500*（t-5）^2+4199/4000
f0=
0.9024
fd0=
-0.3724
2．函数逼近
3．
1.多项式曲线拟合
实例：
用二次多项式拟合下面的数据点
x
1
2
3
y
2
5
10
functionA=multifit（X,Y,m）
N=length（X）;
M=length（Y）;
if（N~=M）
disp（'数据点坐标不匹配!
'）;
return;
end
c（1:
（2*m+1））=0;
b（1:
（m+1））=0;
forj=1:
（2*m+1）
fork=1:
N
c（j）=c（j）+X（k）^（j-1）;
if（j<（m+2））
b（j）=b（j）+Y（k）*X（k）^（j-1）;
end
end
end
C（1,:
）=c（1:
（m+1））;
fors=2:
（m+1）
C（s,:
）=c（s:
（m+s））;
end
A=b'\C;
End
x=1:
3;
y=[2510];
A=multifit（x,y,2）
运行结果
A=
0.12820.32350.8718
2.线性最小二乘拟合
x
1
2
3
4
5
y
1.5
1.8
4
3.4
5.7
矩阵范数
实例：
A=[138;-529;014];
n1=norm（A）
n2=norm（A,1）%列范数
n3=norm（A,inf）%行范数
n4=norm（A,'fro'）%Frobenius范数
运行结果：
n1=
13.5336
n2=
21
n3=
16
n4=
14.1774
3.矩阵条件数
实例：
A=[138;-529;014];
c1=cond（A）
c2=cond（A,1）%1-范数的条件数
c3=cond（A,inf）%行范数的条件数
c4=cond（A,'fro'）%Frobenius范数的条件数
c5=rcond（A）%判断矩阵的病态程度
运行结果：
c1=
40.5924
c2=
85.1053
c3=
61.4737
c4=
42.6696
c5=
0.0118
三．数值微分
1.中点公式法
实例：
中点公式法求一阶导数应用实例，采用中点公式法求函数
的导数。
functiondf=MidPoint（func,x0,h）
ifnargin==2
h=0.1;
elseif（nargin==3&&h==0.0）
disp（'h不能为0!
'）;
return
end
end
y1=subs（sym（func）,findsym（sym（func））,x0+h）;
y2=subs（sym（func）,findsym（sym（func））,x0-h）;
df=（y1-y2）/（2*h）;
df=MidPoint（'sqrt（x）',4）
结果
df=0.2500
2.三点公式法和五点公式法
实例：
采用五点公式法求函数
。
functiondf=ThreePoint（func,x0,type,h）
ifnargin==3
h=0.1;
elseif（nargin==4&&h==0.0）
disp（'h不能为0!
'）
return
end
end
fori=1:
5
y（i）=subs（sym（func）,findsym（sym（func））,x0-2*h+i*h-h）;
end
hd=1/2/h;
switchtype
case1,
df=（-3*y（3）+4*y（4）-y（5））*hd;
case2,
df=（3*y（3）-4*y
（2）+3*y
（1））*hd;
case3,
df=（y（4）-y
（2））*hd;
end
程序
df1=ThreePoint（'sin（x）',2,1）;
df2=ThreePoint（'sin（x）',2,2）;
df3=ThreePoint（'sin（x）',2,3）;
结果
df1=-0.4178
df2=-0.4172
df3=-0.4155
五次
functiondf=FivePoint（func,x0,type,h）
ifnargin==3
h=0.1;
elseif（nargin==4&&h==0.0）
disp（'h不能为0!
'）;
return;
end
end
fori=1:
9
y（i）=subs（sym（func）,findsym（sym（func））,x0-4*h+i*h-h）;
end
hd=1/12/h;
switchtype
case1,
df=（-25*y（5）+48*y（6）-36*y（7）+16*y（8）-3*y（9））*hd;
case2,
df=（-3*y（4）-10*y（5）+18*y（6）-16*y（7）+y（8））*hd;
case3,
df=（y（3）-8*y（4）+8*y（6）-y（7））*hd;
case4,
df=（3*y
（2）+10*y（3）-18*y（4）+16*y（5）-y（6））*hd;
case5,
df=（25*y（5）-48*y（4）+36*y（3）-16*y
（2）+3*y
（1））*hd;
end
程序
df1=FivePoint（'sin（x）',2,1）;
df2=FivePoint（'sin（x）',2,2）;
df3=FivePoint（'sin（x）',2,3）;
df4=FivePoint（'sin（x）',2,4）;
df5=FivePoint（'sin（x）',2,5）;
结果
df1=-0.4161
df2=-0.4161
df3=-0.4161
df4=-0.4161
df5=-0.4161
四．数值积分
1.计算定积分
辛普森法：
function[q,step]=IntSimpon（f,a,b,type,eps）
%被积函数：
f
%积分区间左端点：
a
%积分区间右端点：
b
%辛普森公式类型：
type；
%积分的子区间数：
step
if（type==3&&nargin==4）
disp（'缺少参数'）
end
q=0;
switchtype
case1;
q=（（b-a）/6）*（subs（sym（f）,findsym（sym（f））,a）+...
4*subs（sym（f）,findsym（sym（f））,（a+b）/2+...
subs（sym（f）,findsym（sym（f））,b）））;
step=1;
case2;
q=（（b-a）/8）*（subs（sym（f）,findsym（sym（f）,a）+...
3*subs（sym（f）,findsym（sym（f））,（2*a+b）/3）+...
3*subs（sym（f）,findsym（sym（f））,（a+2*b）/3）+subs（sym（f）,findsym（sym（f））,b）））;
step=1;
case3
n=2;
h=（b-a）/2;
q1=0;
q2=（subs（sym（f）,findsym（sym（f））,a）+subs（sym（f）,findsym（sym（f））,b））/h;
whileabs（q2-q1）>eps
n=n+1;
h=（b-a）/n;
q1=q2;
q2=0;
fori=0:
n-1
x=a+h*i;
x1=x+h;
q2=q2+（h/6）*（subs（sym（f）,findsym（sym（f））,x）+...
4*subs（sym（f）,findsym（sym（f））,（x+x1）/2+...
subs（sym（f）,findsym（sym（f））,x1）））;
end
end
q=q2;
step=n;
end
end
运行结果：
[q,s]=IntSimpon（'sin（x）',0,10,1）
[q,s]=IntSimpon（'sin（x）',0,10,2）
[q,s]=IntSimpon（'sin（x）',0,10,3,0.5）
q=
-6.4487
s=
1
q=
0.7326
s=
1
q=
-0.1088
s=
2
>>[q,s]=IntSimpon（'sin（x）',0,10,3,0.05）
q=
2.0787
s=
19
>>[q,s]=IntSimpon（'sin（x）',0,10,3,0.005）
q=
1.5005
s=
58
五．解线性方程组的方法
1.三对角方程组追赶法
functionx=followup（A,b）
%采用追赶法求线性方程组Ax=b的解
%线性方程组的线性矩阵：
A
%线性方程组中的常数向量：
b
%线性方程组的解:
x
n=rank（A）;
for（i=1:
n）
if（A（i,i）==0）
disp（'Error:
对角有元素为0！
'）;
return;
end
end;
d=ones（n,1）;
a=ones（n-1,1）;
c=ones（n-1）;
for（i=1:
n-1）
a（i,1）=A（i+1,i）;
c（i,1）=A（i,i+1）;
d（i,1）=A（i,i）;
end
d（n,1）=A（n,n）;
%求解Ly=b的解y，解保存在b中，
for（i=2:
n）
d（i,1）=d（i,1）-（a（i-1,1）/d（i-1,1））*c（i-1,1）;
b（i,1）=b（i,1）-（a（i-1,1）/d（i-1,1））*b（i-1,1）;
end
%求解Ux=y的解x,
x（n,1）=b（n,1）/d（n,1）;
for（i=（n-1）:
-1:
1）
x（i,1）=（b（i,1）-c（i,1）*x（i+1,1））/d（i,1）;
end
A=[1-3000;82000;06370;001249;000-45];
b=[1;1;1;1;1];
x=followup（A,b）
x=
0.1923
-0.2692
-0.6923
0.6703
0.7363
2.QR分解法
function[x,Q,R]=qrxq（A,b）
%QR分解法求线性方程组Ax=b的解
%线性方程组的系数矩阵：
A；
%线性方程组的常数矩阵：
b;
%线性方程组的解：
想；
%QR分解后的Q矩阵：
Q；
%QR分解之后的R矩阵：
R；
N=size（A）;
n=N
（1）;
B=A;
B=A;%保存系数矩阵；
A（1:
n,1）=A（1:
n,1）/norm（A（1:
n,i））;%将A的第一列正规化
fori=2:
n
forj=1:
（i-1）
A（1:
n,i）=A（1:
n,i）-dot（A（1:
n,i）,A（1:
n,j））*A（1:
n,j）;
%使A的第i列与前面的所有列正交
end
A（1:
n,i）=A（1:
n,i）/norm（A（1:
n,i））;
%将A的第i列正规化；
end
Q=A;
R=transpose（Q）*B;
x=SolveUpTriangle（R,transpose（Q）*b）;
end
functionx=solveUPTriangle（A,b）
N=size（A）;
n=N
（1）;
dfori=n:
-1:
1
if（is=A（i,（i+1）:
n）*x（（i+1）:
n,1）;
else
s=0;
end
x（i,1）=（b（i）-s）/A（i,i）;
end
A=[832;491;2610];
b=[1;1;1];
[x,Q,R]=qrxq（A,b）
A=[832;491;2610];
b=[1;1;1];
[x,Q,R]=qrxq（A,b）
x=
0.0895
0.0667
0.0421
Q=
0.8729-0.48110.0816
0.43640.6949-0.5715
0.21820.53450.8165
R=
9.16527.85584.3644
-0.00008.01785.0780
0.0000-0.00007.7567
>>
3.超松弛迭代法
function[x,n]=SOR（A,b,x0,w,eps,M）
%采用超松弛迭代法求线性方程组的Ax=b的解
%线性方程组的系数的矩阵：
A;
%线性方程组中的常数矩阵：
吧；
%迭代初始向量；X0;
%松弛因子：
w
%解的精度控制：
eps;
%迭代步数控制：
M
%线性方程组的解：
想；
%求出所需精度的解的实际迭代次数；n
ifnargin==4
eps=1.0e-6;
M=10000;
elseifnargin==5;
M=10000;
end
if（w<=0||w>2）
error
return;
end
D=diag（diag（A））;
L=-tril（A,-1）;
U=-triu（A,1）;
B=inv（D-L*w）*（（1-w）*D+w*U）;
f=w*inv（（D-L*w））*b;
x=x0;
n=0;
tol=1;
whiletol>=eps
x=B*x0+f;
n=n+1;
tol=norm（x-x0）;
x0=x;
if（n>=M）
diap（'Warning:
迭代次数太多，可能不收敛！
'）;
return;
end
end
end
运行结果：
A=[0.680.010.12;0.03-0.54-0.05;0.20.080.74];
b=[1;1;1];
x0=[0;0;0];
[x,n]=SOR（A,b,x0,1.07）
x=
1.2851
-1.8924
1.2086
n=
7
4.雅可比迭代法
函数调用
function[x,n]=jacobi（A,b,x0,eps,M）
%采用雅可比迭代法求解线性方程组Ax=b的解
%线性方程组的系数矩阵:
A
%线性方程组的常数向量:
b
%迭代初始向量:
x0
%解的精度控制:
eps
%迭代步数控制:
M
%线性方程组的解；x
%求出所需精度的解实际的迭代步数：
n
ifnargin==3
eps=1.0e-6;
M=10000;
elseifnargin==4
M=10000;
end
D=diag（diag（A））;%求A的对角矩阵
L=-tril（A,-1）;%求A的下三角阵
U=-triu（A,1）;%求A的上三角阵
B=D\（L+U）;
f=D\b;
x=x0;
n=0;%迭代次数
tol=1;%前后两次迭代结果误差
whiletol>=eps
x=B*x0+f;%迭代公式
n=n+1;
tol=norm（x-x0）;
x0=x;
if（n>=M）
dosp（'Waring:
迭代次数太多，可能不收敛!
'）;
return;
end
end
运行程序；
A=[0.98-0.05-0.02;-0.04-0.90.07;-0.020.090.94];
b=[1;1;1];
x0=[0;0;0];
[x,n]=jacobi（A,b,x0）
结果：
x=
0.9904
-1.0628
1.1867
n=
8
5.高斯-赛德尔迭代法
函数调用
function[x,n]=gauseidel（A,b,x0,eps,M）
%采用高斯-赛德尔迭代法求解线性方程组Ax=b的解
%线性方程组的系数矩阵:
A
%线性方程组的常数向量:
b
%迭代初始向量:
x0
%解的精度控制:
eps
%迭代步数控制:
M
%线性方程组的解；x
%求出所需精度的解实际的迭代步数：
n
ifnargin==3
eps=1.0e-6;
M=10000;
elseifnargin==4
M=10000;
end
D=diag（diag（A））;%求A的对角矩阵
L=-tril（A,-1）;%求A的下三角阵
U=-triu（A,1）;%求A的上三角阵
G=（D-L）\U;
f=（D-L）\b;
x=x0;
n=0;%迭代次数
tol=1;%前后两次迭代结果误差
whiletol>=eps
x=G*x0+f;%迭代公式
n=n+1;
tol=norm（x-x0）;
x0=x;
if（n>=M）
dosp（'Waring:
迭代次数太多，可能不收敛!
'）;
return;
end
end
程序
A=[0.98-0.05-0.02;-0.04-0.90.07;-0.020.090.94];
b=[1;1;1];
x0=[0;0;0];
[x,n]=gauseidel（A,b,x0）
结果
x=
0.9904
-1.0628
1.1867
n=
5
6.牛顿插值法
函数
function[f,f0]=Newton（x,y,x0）
symst;
if（length（x）==length（y））
n=length（x）;
c（1:
n）=0.0;
else
disp（'xºÍyµÄάÊý²»ÏàµÈ£¡'）;
return;
end
f=y
（1）;
y1=0;
l=1;
for（i=1:
n-1）
升级会员 