Solving equation in a for loop with data
2 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Jonas Damsbo
el 6 de En. de 2020
Comentada: Jonas Damsbo
el 7 de En. de 2020
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
0 comentarios
Respuesta aceptada
David Hill
el 7 de En. de 2020
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 comentarios
Más respuestas (1)
Walter Roberson
el 7 de En. de 2020
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.
0 comentarios
Ver también
Categorías
Más información sobre Function Handles en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!