Solving equation in a for loop with data

1 次查看(过去 30 天)
I have this code with a for loop where I want to determine h(i).
Is there any way I can solve this equation for h?
h is the ice thikness. x is distance from a glacier edge, M is the glacier sliced in 100 elements. y is gravity data and the equation in the for loop is the equation for gravity anomalies.
%Data
x = [535 749 963 1177 1391 1605 1819 2033 2247 2461 2675 2889]; %Distance from edge (m)
y = [-15.0 -24.0 -31.2 -36.8 -40.8 -42.7 -42.4 -40.9 -37.3 -31.5 -21.8 -12.8].*10.^-5; %Anomaly (mgal)
M = 1:34.2:3420; %M regular elements
dx = 34.2; %M is dx wide
%Plot data
figure(1);
plot(x,y,'*');
xlabel('Distance from west edge (m)');
ylabel('Bouguer anomaly (mgal)');
%The inverse problem
%(1) Perform discretization
G = 6.67*10^(-11); %Gravity constant in [m3/ kg*s^2]
del_rho = -1700; %Density contrast
d = 0.1;
for j = 1:12 %We have 12 measurements
for i = 1:100 %100 M regular elements
y(j) = G.*del_rho.*ln((M(i)-x(j).^2 + h(i).^2)./(M(i) - x(j)).^2 + d).*dx
end
end

采纳的回答

David Hill
David Hill 2020-1-7
x = [535 749 963 1177 1391 1605 1819 2033 2247 2461 2675 2889]; %Distance from edge (m)
y = [-15.0 -24.0 -31.2 -36.8 -40.8 -42.7 -42.4 -40.9 -37.3 -31.5 -21.8 -12.8].*10.^-5; %Anomaly (mgal)
M = 1:34.2:3420; %M regular elements
dx = 34.2; %M is dx wide
%Plot data
figure(1);
plot(x,y,'*');
xlabel('Distance from west edge (m)');
ylabel('Bouguer anomaly (mgal)');
%The inverse problem
%(1) Perform discretization
G = 6.67*10^(-11); %Gravity constant in [m3/ kg*s^2]
del_rho = -1700; %Density contrast
d = 0.1;
h=zeros(1,1200);
for j = 1:12 %We have 12 measurements
for i = 1:100 %100 M regular elements
h((j-1)*100+i)=sqrt(exp(y(j)/(G*del_rho*dx))*(M(i)-x(j))^2 + x(j)^2 - M(i));%you can directly solve for h but there are 1200 different values
end
end
semilogy(1:1200,h);%there are 100 values of h for each of the 12 measurements, you could reshape(h,12,100) to get rows for each measurement.
I am not sure this is what you want, but you can easily solve for h in your equation.
  2 个评论
Jonas Damsbo
Jonas Damsbo 2020-1-7
Maybe Im not sure. My problem is an inverse problem. I have the 12 measurements g(j) and these 12 measurements are located at x (j = 1:12) on the x-axis. The equation for g in my for loop is then M sums (the log part), where each sum has the height h(i). Here x (i) is a random selected value out of the x-axis - e.g. 100 step. For each sum, there is then an estimated height h (i). And this is the one I want to get out as a vector. Not sure if this is possible with an algorithm? Or whether g (j) should be the same size as h?
Jonas Damsbo
Jonas Damsbo 2020-1-7
That the equation where I want the estimates of h.

请先登录,再进行评论。

更多回答(1 个)

Walter Roberson
Walter Roberson 2020-1-7
syms h
syms d dx G M positive
syms del_rho real
assume(del_rho < 0)
syms x y real
eqn = y == G.*del_rho.*log((M-x.^2 + h.^2)./(M - x).^2 + d).*dx;
sol = solve(eqn,h);
There are two solutions, which are +/-
(M^2*exp(y/(G*del_rho*dx)) - M^2*d - M - d*x^2 + x^2*exp(y/(G*del_rho*dx)) + x^2 + 2*M*d*x - 2*M*x*exp(y/(G*del_rho*dx)))^(1/2)
You can subs() in the numeric arrays to get the vector of solutions.

类别

Help CenterFile Exchange 中查找有关 Oceanography and Hydrology 的更多信息

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by