how to solve system of equation

5 visualizaciones (últimos 30 días)
Priyank Pathak
Priyank Pathak el 11 de Sept. de 2021
Comentada: Priyank Pathak el 12 de Sept. de 2021
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% I have two above given equations in the picture.
% i used these the given code, but it is taking too much time to run and showing some error in "solve". please suggest, what i can do. Or if possible , please write some new code.
%%%%%%%%%%%%%%%%%%%%%%%%%%%
E=xlsread('2131toposhg','c1:c90601');
N=xlsread('geoid anomaly latlong','c1:c90601');
L0=2400; %load in metres
r_c=2670; % density of crust in kg/m3
r_w=1030; % density of water in kg/m3
r_a=3200; % density of asthenosphere in kg/m3
r_m=3250; % density of lithosphere mantle in kg/m3
Zmax=300000; %compensation level, zmax in metres
g=9.8; % gravity in m/s2
r_ac=r_a-r_c; % density contrast b/w asthenosphere and crust
r_cw=r_c-r_w; % density contrast b/w crust and water
r_wc=r_w-r_c; % density contrast b/w water and crust
r_ma=r_m-r_a; % density contrast b/w mantle lithosphere and asthenosphere
r_mc=r_m-r_c; % density contrast b/w mantle lithosphere and crust
ur_ac=1/(r_ac);
ZL= 130000; %in metres
ZC= 47000; % in metres
K=r_a*L0+E*r_cw;
C2=3.14*6.67*10^(-11);
E2=E.^2;
ZC2=ZC.^2;
ZL2=ZL.^2;
Zmax2=Zmax^2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N0=-(C2*(E2*r_wc+ (ZC2-E2)*r_c+(ZL2-ZC2)*r_m+ (Zmax2-ZL2)*r_a))./ 9.8 -N ;
N22=-(N0+N)*9.8/C2;
syms zL zC
eqn1= zC-(K+zL*r_ma)./(r_mc)==0;
eqn2= C2*(E2*r_wc+ (zL.^2-E2)*r_c+(zL.^2-zC.^2)*r_m+ (Zmax2-zL.^2)*r_a)==N22 ;
sol = solve([eqn1, eqn2], [zL, zC]);
zLSol = sol.zL
zCSol = sol.zC

Respuestas (1)

Walter Roberson
Walter Roberson el 11 de Sept. de 2021
E=xlsread('2131toposhg','c1:c90601');
N=xlsread('geoid anomaly latlong','c1:c90601');
vectors
K=r_a*L0+E*r_cw;
vector because E is a vector
E2=E.^2;
vector becaue E is a vector
N0=-(C2*(E2*r_wc+ (ZC2-E2)*r_c+(ZL2-ZC2)*r_m+ (Zmax2-ZL2)*r_a))./ 9.8 -N ;
vector because E2 is a vector (because E is a vector) and because N is a vector. Subtraction should be valid because E2 and N should be the same size.
N22=-(N0+N)*9.8/C2;
vector because N0 is a vector (becaue E2 is a vector because E is a vector)
syms zL zC
zL and zC are scalars.
eqn1= zC-(K+zL*r_ma)./(r_mc)==0;
vector because K is a vector (becaue E is a vector). Involves the scalar symbolic variables zC and zL
eqn2= C2*(E2*r_wc+ (zL.^2-E2)*r_c+(zL.^2-zC.^2)*r_m+ (Zmax2-zL.^2)*r_a)==N22 ;
vector on the left side because E2 is a vector (because E is a vector) . Vector on the right hand side because N22 is a vector (because N0 is a vector because E2 is a vector because E is a vector). The left and right side of the == should be the same size. Also, eqn2 should be the same size as eqn1 . Involves the scalar symbolic variables zL and zC
sol = solve([eqn1, eqn2], [zL, zC]);
eqn1 and eqn2 are vector the same length, both involving the scalar symbolic variables zL and zC .
You have 90601 + 90601 equations created by the [eqn1, eqn2], all of which involve the scalar symbolic variables zC and zL. You are asking to solve that system of 181202 equations for the two scalar variables zL and zC . That is a rather overdetermined system. It is very unlikely that you will be able to find a single zL and single zC that will solve a 181202 equations at the same time.
solve() is for solving systems of simultaneous equations to find the values of the variables that make all of the equations true simultaneously.
solve() is not vectorized -- if you give it multiple rows or multiple columns of equations, it will not attempt to solve each row or column independently.
You need to take one of two approaches:
  1. Make zC and zL symbolic vectors of size 90601 * 1 each, so that effectively you would be solving 181202 equations for 181202 variables. If you do this then the solver output, sol would be a struct with 181202 fields. Probably not optimal. OR
  2. arrayfun(@(E1, E2) solve([E1, E2], [zL, zC]), eqn1, eqn2, 'uniform', 0) . This would return a cell array with 90601 entries in which each entry was a struct of solutions.
  5 comentarios
Walter Roberson
Walter Roberson el 12 de Sept. de 2021
In the below, you might need to add the option 'ReadVariableNames', false to the readtable() calls.
If your version of MATLAB is too old to have readtable() then you will need to switch back to xlsread() -- but xlsread will definitely not be able use the URLs for the names of the files to read.
top_file = 'https://www.mathworks.com/matlabcentral/answers/uploaded_files/735774/2131toposhg.xlsx';
geoid_file = 'https://www.mathworks.com/matlabcentral/answers/uploaded_files/735784/geoid%20anomaly%20latlong.xlsx';
syms N E real
L0=2400; %load in metres
r_c=2670; % density of crust in kg/m3
r_w=1030; % density of water in kg/m3
r_a=3200; % density of asthenosphere in kg/m3
r_m=3250; % density of lithosphere mantle in kg/m3
Zmax=300000; %compensation level, zmax in metres
g=9.8; % gravity in m/s2
r_ac=r_a-r_c; % density contrast b/w asthenosphere and crust
r_cw=r_c-r_w; % density contrast b/w crust and water
r_wc=r_w-r_c; % density contrast b/w water and crust
r_ma=r_m-r_a; % density contrast b/w mantle lithosphere and asthenosphere
r_mc=r_m-r_c; % density contrast b/w mantle lithosphere and crust
ur_ac=1/(r_ac);
ZL= 130000; %in metres
ZC= 47000; % in metres
K=r_a*L0+E*r_cw;
C2=3.14*6.67*10^(-11);
E2=E.^2;
ZC2=ZC.^2;
ZL2=ZL.^2;
Zmax2=Zmax^2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N0=-(C2*(E2*r_wc+ (ZC2-E2)*r_c+(ZL2-ZC2)*r_m+ (Zmax2-ZL2)*r_a))./ 9.8 -N ;
N22=-(N0+N)*9.8/C2;
syms zL zC positive
eqn1= zC-(K+zL*r_ma)./(r_mc)==0;
eqn2= C2*(E2*r_wc+ (zL.^2-E2)*r_c+(zL.^2-zC.^2)*r_m+ (Zmax2-zL.^2)*r_a)==N22 ;
sol = solve([eqn1, eqn2], [zL, zC], 'returnconditions', true)
sol = struct with fields:
zL: [2×1 sym] zC: [2×1 sym] parameters: [1×0 sym] conditions: [2×1 sym]
simplify(sol.conditions)
ans = 
E_actual = readtable(top_file, 'range', 'c1:c90601');
N_actual = readtable(geoid_file, 'range', 'c1:c90601');
HowManyForDemo = 20;
sol_actual = subs(sol, {E, N}, {E_actual{1:HowManyForDemo, 1}.', N_actual{1:HowManyForDemo,1}.'})
sol_actual = struct with fields:
zL: [2×20 sym] zC: [2×20 sym] parameters: [1×0 sym] conditions: [2×20 sym]
vpa(sol_actual.zL(2,:))
ans = 
vpa(sol_actual.zC(2,:))
ans = 
Priyank Pathak
Priyank Pathak el 12 de Sept. de 2021

Iniciar sesión para comentar.

Etiquetas

Community Treasure Hunt

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

Start Hunting!

Translated by