物探数据反演上机实验.docx
《物探数据反演上机实验.docx》由会员分享,可在线阅读,更多相关《物探数据反演上机实验.docx(27页珍藏版)》请在冰豆网上搜索。
![物探数据反演上机实验.docx](https://file1.bdocx.com/fileroot1/2023-2/26/767ce524-fa7c-478e-a63d-ad44c218c454/767ce524-fa7c-478e-a63d-ad44c218c4541.gif)
物探数据反演上机实验
中南大学实验报告
物探数据处理与反演上机实验
班级:
姓名:
学号:
指导老师:
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)