function [alpha,N,U]=youxianchafen2(r1,r2,up,under,num,deta)
%[alpha,N,U]=youxianchafen2(a,r1,r2,up,under,num,deta)
%該函數(shù)用有限差分法求解有兩種介質(zhì)的正方形區(qū)域的二維拉普拉斯方程的數(shù)值解
%函數(shù)返回迭代因子、迭代次數(shù)以及迭代完成后所求區(qū)域內(nèi)網(wǎng)格節(jié)點處的值
%a為正方形求解區(qū)域的邊長
%r1,r2分別表示兩種介質(zhì)的電導(dǎo)率
%up,under分別為上下邊界值
%num表示將區(qū)域每邊的網(wǎng)格剖分個數(shù)
%deta為迭代過程中所允許的相對誤差限
n=num+1; %每邊節(jié)點數(shù)
U(n,n)=0; %節(jié)點處數(shù)值矩陣
N=0; %迭代次數(shù)初值
alpha=2/(1+sin(pi/num));%超松弛迭代因子
k=r1/r2; %兩介質(zhì)電導(dǎo)率之比
U(1,1:n)=up; %求解區(qū)域上邊界第一類邊界條件
U(n,1:n)=under; %求解區(qū)域下邊界第一類邊界條件
U(2:num,1)=0;U(2:num,n)=0;
for i=2:num
U(i,2:num)=up-(up-under)/num*(i-1);%采用線性賦值對上下邊界之間的節(jié)點賦迭代初值
end
G=1;
while G>0 %迭代條件：不滿足相對誤差限要求的節(jié)點數(shù)目G不為零
Un=U; %完成第n次迭代后所有節(jié)點處的值
G=0; %每完成一次迭代將不滿足相對誤差限要求的節(jié)點數(shù)目歸零
for j=1:n
for i=2:num
U1=U(i,j); %第n次迭代時網(wǎng)格節(jié)點處的值
if j==1 %第n+1次迭代左邊界第二類邊界條件
U(i,j)=1/4*(2*U(i,j+1)+U(i-1,j)+U(i+1,j));
end
if (j>1)&&(j U2=1/4*(U(i,j+1)+ U(i-1,j)+ U(i,j-1)+ U(i+1,j));
U(i,j)=U1+alpha*(U2-U1); %引入超松弛迭代因子后的網(wǎng)格節(jié)點處的值
end
if i==n+1-j %第n+1次迭代兩介質(zhì)分界面（與網(wǎng)格對角線重合）第二類邊界條件
U(i,j)=1/4*(2/(1+k)*(U(i,j+1)+U(i+1,j))+2*k/(1+k)*(U(i-1,j)+U(i,j-1)));
end
if j==n %第n+1次迭代右邊界第二類邊界條件
U(i,n)=1/4*(2*U(i,j-1)+U(i-1,j)+U(i+1,j));
end
end
end
N=N+1 %顯示迭代次數(shù)
Un1=U; %完成第n+1次迭代后所有節(jié)點處的值
err=abs((Un1-Un)./Un1);%第n+1次迭代與第n次迭代所有節(jié)點值的相對誤差
err(1,1:n)=0; %上邊界節(jié)點相對誤差置零
err(n,1:n)=0; %下邊界節(jié)點相對誤差置零
G=sum(sum(err>deta))%顯示每次迭代后不滿足相對誤差限要求的節(jié)點數(shù)目G
end
標簽：
有限差分
上傳時間：
2018-07-13
上傳用戶：Kemin
%球體
close all;
G=6.67e-11;
R=2;%球體半徑
p=4.0;%密度
D=10.0;%深度
M=(4/3)*pi*R^3*p;%質(zhì)量
x=-20:1:20;
g=G*M*D./((x.^2+D^2).^(3/2));
Vxz=-3*G*M*D.*x./((x.^2+D^2).^(5/2));
Vzz=G*M.*(2*D^2-x.^2)./((x.^2+D^2).^(5/2));
Vzzz=3*G*M.*(2*D^2-3.*x.^2)./((x.^2+D^2).^(7/2));
subplot(2,2,1)
plot(x,g,'k-');
xlabel('水平距離(m)');
ylabel('重力異常值');
title('球體重力異常Δg');
grid on
subplot(2,2,2)
plot(x,Vxz);
xlabel('水平距離(m)');
ylabel('導(dǎo)數(shù)值');
title('Vxz');
grid on
subplot(2,2,3)
plot(x,Vzz);
xlabel('水平距離(m)');
ylabel('導(dǎo)數(shù)值');
title('Vzz');
grid on
subplot(2,2,4);
plot(x,Vzzz);
xlabel('水平距離(m)');
ylabel('導(dǎo)數(shù)值');
title('Vzzz');
grid on
%%
%水平圓柱體
close all
G=6.67e-11;
p=10.0;%線密度
D=100.0;%深度
x=-200:1:200;
g=G*2*p*D./(x.^2+D^2);
Vxz=4*G*p*D.*x./(x.^2+D^2).^2;
Vzz=2*G*p.*(D^2-x.^2)./(x.^2+D^2).^2;
Vzzz=4*G*p.*(D^2-3.*x.^2)./((x.^2+D^2).^3);
subplot(2,2,1)
plot(x,g,'k-');
xlabel('水平距離(m)');
ylabel('重力異常值');
title('水平圓柱體重力異常Δg');
grid on
subplot(2,2,2)
plot(x,Vxz);
xlabel('水平距離(m)');
ylabel('導(dǎo)數(shù)值');
title('Vxz');
grid on
subplot(2,2,3)
plot(x,Vzz);
xlabel('水平距離(m)');
ylabel('導(dǎo)數(shù)值');
title('Vzz');
grid on
subplot(2,2,4);
plot(x,Vzzz);
xlabel('水平距離(m)');
ylabel('導(dǎo)數(shù)值');
title('Vzzz');
grid on
%%
%垂直臺階
G=6.67e-11;
p=4.0;%密度
h1=50.0;%下層深度
h2=40.0;%上層深度
x=-100:1:100;
g=G*p.*(pi*(h1-h2)+x.*log((x.^2+h1^2)./(x.^2+h2^2))+2*h1.*atan(x./h1)-2*h2.*atan(x./h2));
Vxz=G*p.*log((h1^2+x.^2)./(h2^2+x.^2));
Vzz=2*G*p.*atan((x.*(h1-h2))./(x.^2+h1*h2));
Vzzz=2*G*p.*x*(h1^2-h2^2)./((h1^2+x.^2).*(x.^2+h2^2));
subplot(2,2,1)
plot(x,g,'k-');
xlabel('水平距離(m)');
ylabel('重力異常值');
title('垂直臺階重力異常Δg');
grid on
subplot(2,2,2)
plot(x,Vxz);
xlabel('水平距離(m)');
ylabel('導(dǎo)數(shù)值');
title('Vxz');
grid on
subplot(2,2,3)
plot(x,Vzz);
xlabel('水平距離(m)');
ylabel('導(dǎo)數(shù)值');
title('Vzz');
grid on
subplot(2,2,4);
plot(x,Vzzz);
xlabel('水平距離(m)');
ylabel('導(dǎo)數(shù)值');
title('Vzzz');
grid on
%%
%傾斜臺階
G=6.67e-11;
p=4.0;%密度
h1=50.0;%下層深度
h2=40.0;%上層深度
a=pi/6;%傾斜角度
x=-500:1:500;
g=G*p.*(pi*(h1-h2)+2*h1.*atan((x+h1*cot(a))./h1)-2*h2.*atan((x+h2*cot(a))./h1)+x.*sin(a)^2.*log(((h1+x.*sin(a).*cos(a)).^2+x.^2.*sin(a)^4)./((h2+x.*(sin(a)*cos(a))).^2+x.^2.*sin(a)^4)));
Vxz=G*p.*(sin(a)^2.*log(((h1*cot(a)+x).^2+h1^2)./((h2*cot(a)+x).^2+h2^2))-2*sin(2*a).*(atan((h1/sin(a)+x.*cos(a))./(x.*sin(a)))-atan((h2/sin(a)+x.^cos(a))./(sin(a).*x))));
Vzz=G*p.*(0.5*sin(2*a)^2.*log(((h1*cot(a)+x).^2+h1^2)./((h2*cot(a)+x).^2+h2^2))+2*sin(a)^2.*(atan((h1/sin(a)+x.*cos(a))./(x.*sin(a)))-atan((h2/sin(a)+x.*cos(a))./(x.*sin(a)))));
Vzzz=2*G*p*sin(a)^2.*((x+2*h2*cot(a))./((h2*cot(a)+x).^2+h2^2)-(x+2*h1*cot(a))./((h1*cot(a)+x).^2+h1^2));
subplot(2,2,1)
plot(x,g,'k-');
xlabel('水平距離(m)');
ylabel('重力異常值');
title('傾斜臺階重力異常Δg');
grid on
subplot(2,2,2)
plot(x,Vxz);
xlabel('水平距離(m)');
ylabel('導(dǎo)數(shù)值');
title('Vxz');
grid on
subplot(2,2,3)
plot(x,Vzz);
xlabel('水平距離(m)');
ylabel('導(dǎo)數(shù)值');
title('Vzz');
grid on
subplot(2,2,4);
plot(x,Vzzz);
xlabel('水平距離(m)');
ylabel('導(dǎo)數(shù)值');
title('Vzzz');
grid on
%%
%鉛錘柱體
G=6.67e-11;
p=4.0;%密度
h1=50.0;%下層深度
h2=40.0;%上層深度
a=3;%半徑
x=-500:1:500;
g=G*p.*((x+a).*log(((x+a).^2+h1^2)./((x+a).^2+h2^2))-(x-a).*log(((x-a).^2+h1^2)./((x-a).^2+h2^2))+2*h1.*(atan((x+a)./h1)-atan((x-a)./h1))-2*h2.*(atan((x+a)./h2)-atan((x-a)./h2)));
Vxz=G*p.*log((((x+a).^2+h1^2).*((x-a).^2+h2^2))./(((x+a).^2+h2^2).*((x-a).^2+h1^2)));
Vzz=2*G*p.*(atan(h1./(x+a))-atan(h2./(x+a))-atan(h1./(x-a))+atan(h2./(x-a)));
Vzzz=2*G*p.*((x+a)./((x+a).^2+h2^2)-(x+a)./((x+a).^2+h1^2)-(x-a)./((x-a).^2+h2^2)+(x-a)./((x-a).^2+h1^2));
subplot(2,2,1)
plot(x,g,'k-');
xlabel('水平距離/m')
ylabel('重力異常值')
title('鉛垂柱體重力異常')
grid on
subplot(2,2,2)
plot(x,Vxz);
xlabel('水平距離(m)');
ylabel('導(dǎo)數(shù)值');
title('Vxz');
grid on
subplot(2,2,3)
plot(x,Vzz);
xlabel('水平距離(m)');
ylabel('導(dǎo)數(shù)值');
title('Vzz');
grid on
subplot(2,2,4);
plot(x,Vzzz);
xlabel('水平距離(m)');
ylabel('導(dǎo)數(shù)值');
title('Vzzz');
grid on
標簽：
MATLAB
重力
正
程序
上傳時間：
2019-05-10
上傳用戶：xiajiang
