- 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
- 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.

# how to solve system of equation

5 views (last 30 days)

Show older comments

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 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

##### 0 Comments

### Answers (1)

Walter Roberson
on 11 Sep 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:

##### 5 Comments

Walter Roberson
on 12 Sep 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)

simplify(sol.conditions)

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}.'})

vpa(sol_actual.zL(2,:))

vpa(sol_actual.zC(2,:))

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!