Solving equation in a for loop with data

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

Answer 1

[David Hill]

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.

[Jonas Damsbo]

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]

That the equation where I want the estimates of h.

Answer 2

[Walter Roberson]

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.