how to solve system of equation

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

Answer 1

[Walter Roberson]

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.

[Walter Roberson]

Not tested because I do not have the file.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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');

Error using xlsread (line 136)
Unable to open file '2131toposhg'.
File '2131toposhg' not found.

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_cell = arrayfun(@(E1,E2) solve([E1, E22], [zL, zC]), eqn1, eqn2, 'uniform', 0);
zLSol = cellfun(@(C) C.zL, sol_cell, 'uniform', 0)
zCSol = cellfun(@(C) C.zC, sol_cell, 'uniform', 0)

[Walter Roberson]

There is a pair of solutions at each point -- roots of a quadartic. In the values I looked at, one of the pair (the first row) was always negative, and the other was always positive -- not exact negatives of each other, but very similar.

... Watch out for the HowManyForDemo = 20; in the below -- I only had it display a subset for demontration purposes.

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 = readmatrix(top_file, 'range', 'c1:c90601');
N_actual = readmatrix(geoid_file, 'range', 'c1:c90601');
HowManyForDemo = 20;
sol_actual = subs(sol, {E, N}, {E_actual(1:HowManyForDemo).', N_actual(1:HowManyForDemo).'})

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]

@Walter Roberson, When i am running it shows error in "read matrix".

the command window comment.-----------

Undefined function or variable 'readmatrix'.

Error in lbprogram (line 36)

E_actual = readmatrix(top_file, 'range', 'c1:c90601');

[Walter Roberson]

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]

@Walter Roberson Thank you