How do I solve lhs=rhs by iteration?

For my research I want to use decompaction and ultimately backstripping to solve subsidence of geological layers in the subsurface. I could do this for one layer by hand, but I want to write a script that will automatically solve this for all layers in my borehole data.

For the first step in decompaction, I want to solve the following function:

y2'-y1'=y2-y1-phi0/c*{exp-(c*y1)-exp(-c*y2)}+phi0/c*{exp(-c*y1')-exp(-c*y2')

The following variables I already have: y1',y1 and y2. I want to solve for y2', but using test-data (for which I have the answers) I keep getting the wrong answers.

(y1 is the initial top of the layer, y2 is the base and y1' is the shallower depth for decompaction) For simplicity purposes, I use y1'=0.

My test-data are: y1= 3000 y2= 4000 y1'=0 phi0=0.70 c=0.71

the answer for this test-data for my script should be y2'=1620.

can anyone help me solve my problem?

I am not proficient enough in Matlab to solve this problem myself, but I would love to learn.

[Star Strider]

Do the primes (') signify derivatives with respect to time, or matrix or vector complex conjugate transposes?

[Annemijn van den Berg]

No, the primes (') signify that these are new top and base values for the same layer. could also be called y1_1 and y2_1.

[Walter Roberson]

y2'-y1'=y2-y1-phi0/c*{exp-(c*y1)-exp(-c*y2)}+phi0/c*{exp(-c*y1')-exp(-c*y2')

has bad brackets. At the very least it is missing a } . The exp should not be followed by a -

Answer 1

[Star Strider]

Geology is regrettably not an area of my expertise, so I am not familiar with what you are doing.

It might be best to approach this as an optimization problem. Since you have all but one variable, you can use either fzero or fsolve to solve for y2' (that I call ‘y2p’ here).

This is returning undefined values at the initial point, so I ask you to be certain the parentheses are correct, since they are not precise in your equation, and I might have placed them incorrectly in mine:

y1   = 3000;
y2   = 4000;
y1p  = 0;
phi0 = 0.70;
c    = 0.7;

fcn = @(y2p) y2-y1-phi0/(c*(exp(-c*y1)-exp(-c*y2)))+phi0/(c*(exp(-c*y1p)-exp(-c*y2p))) - y2p + y1p;
y2p = fzero(fcn, 2E+3);

I rearranged your equation so that it is implicit (equals zero) for the fzero and fsolve functions, since they require that.

I used fzero here because everyone has it. The Optimization Toolbox fsolve function is more robust. With correctly-placed parentheses, that should produce the correct result.

I will help as I can to get this working.

[Annemijn van den Berg]

Okay, the equation looks good, but then I get the following error:

Error using fzero (line 328)
Function value at starting guess must be finite
and real.

What did I do wrong?

[Star Strider]

That was the same error I got.

The problem is likely that both exp(-c*y1) and exp(-c*y2) are 0 at the outset. That means that the first term phi0/(c*(exp(-c*y1)-exp(-c*y2))) is infinite. The ‘y1’ and ‘y2’ values are numerically too large to produce a non-zero result.

Even letting ‘c’ vary as a parameter results in a solution that does not vary much from the initial estimate. I also tried this with the Symbolic Math Toolbox (that has greater numeric precision). It could not find a solution.

I am willing to entertain any ideas you may have.

[Annemijn van den Berg]

I tried the following, which obviously got me the wrong answer..

Phi0=[0.63; 0.49; 0.70; 0.40; 0.20; 0.05];
c=[0.51; 0.27; 0.71; 0.60; 0.60; 0.20];
y1= 3000;
y2= 4000;
y1_1=zeros(1,5001);
y2_2=[0:5000];
lhs= y2_2-y1_1;
rhs= y2-y1-Phi0(3)./c(3).*(exp(-c(3).*y1)-exp(-c(3).*y2))+Phi0(3)./c(3).*(exp(-c(3).*y1_1)-exp(-c(3).*y2_2));

P=intersect(lhs, rhs);

The intersection here is just the difference between y1 and y2 (so the thickness)

I think modifying what I did here should work, but I just can't seem to find what I did wrong.

Thank you for wanting to help me!

[Star Strider]

My pleasure.

Running this code with your new equation and new constants:

fcn = @(y2_2) y2-y1-Phi0(3)./c(3).*(exp(-c(3).*y1)-exp(-c(3).*y2))+Phi0(3)./c(3).*(exp(-c(3).*y1_1)-exp(-c(3).*y2_2)) - (y2_2-y1_1);
y2_2 = fsolve(fcn, 2000)

y2_2 =
         1001

This equation works.

[Annemijn van den Berg]

Yes the equation seems to work, but for the test-data I know the answer and the answer should be 1612. Maybe Matlab isn't going to help me solve this for my real data...

[Star Strider]

I trust MATLAB’s calculations.

Your equation as written appears to be consistent with respect to the numeric and analytic solutions. If it is not producing the result you want, carefully review your equation and known values to be certain all of them are correct.

I will help you with this if I know what your equation actually is. Without that, I cannot be certain it is coded correctly. If the equation and constants are correct as written here, are you certain that the constants and the 1612 value are also correct?

[Annemijn van den Berg]

I actually found a script online from someone who was working on something similar and found a way to make it work!

Thanks for your help though, it is very much appreciated.

[Star Strider]

My pleasure.

Please post the script you found here. I would like to see it.

[Annemijn van den Berg]

http://www.blackwellpublishing.com/allen/practicalexercises/exercise9_1.m

this script takes care of the decompaction as well as the backstripping in one go.

Answer 2

[Walter Roberson]

If the equation is

y2_new-y1_new = y2-y1-phi0*(exp(-c*y1)-exp(-c*y2))/c+phi0*(exp(-c*y1_new)-exp(-c*y2_new))/c;

then this has a solution f

 y2_new = (LambertW(-phi0*exp(-phi0*exp(-c*y1_new)+phi0*exp(-c*y1)-exp(-c*y2)*phi0+(y1-y1_new-y2)*c))+phi0*exp(-c*y1_new)-phi0*exp(-c*y1)+exp(-c*y2)*phi0+(-y1+y2+y1_new)*c)/c

With real inputs, LambertW often has exactly two real solutions, one denoted as LambertW(-1,x) and the other as LambertW(0,x) which gets abbreviated to LambertW(x) .

With the coefficients you gave in the original Question, the corresponding real solutions are approximately -9.7643 and 1000.9859

[Annemijn van den Berg]

Yes the equation seems to work, but for the test-data I know the answer and the answer should be 1612. Maybe Matlab isn't going to help me solve this for my real data...