Recursive formulas - Numerical method need

Hello to everyone,i'm trying to obtain the pressures under a footing applying the procedure in Complete Analytical Solution for Linear Soil Pressure Distribution under Rigid Rectangular Spread Footings by John Bellos and Nikolaos Bakas, but i'm stuck in this set of equations

the unknown variables are xn and yn, there is any numerical method to solve this kind of problem?

[Torsten]

2 equations for 2 unknowns for each value of n. Use MATLAB's fsolve, e.g.

Fixpoint iteration might also work.

[Roberto Enrique Pinto Villegas]

Thank you, i will try it.

[Walter Roberson]

It would help if you could provide sample values for the known variables, for us to test with.

[Roberto Enrique Pinto Villegas]

I send an email to John Bellos asking how he done the calculations, he say We used recursion to find the solution, starting from certain values of xn and yn, we calculate new xn and yn and so on so i tried the following

    function [xn,yn]=coordinates(lx,ly,ex,ey,rex,rey)
    syms xn yn
    
    eqx=4*lx*(1/2-rex)*(1-(1-lx/xn)^3-(1-ly/yn)^3)/(1-(1-lx/xn)^4-(1-ly/yn)^4-4*...
        (1-lx/xn)^3*lx/xn);
    eqy=4*ly*(1/2-rey)*(1-(1-lx/xn)^3-(1-ly/yn)^3)/(1-(1-lx/xn)^4-(1-ly/yn)^4-4*...
        (1-ly/yn)^3*ly/yn);
    
    xn=lx;
    yn=ly;
    
    x1=subs(eqx);
    y1=subs(eqy);
    
    e1=xn-x1;
    e2=yn-y1;
    
    while ((e1<=-0.01) || (e1>=0.01)) && ((e2<=-0.01) || (e2>=0.01))
        xn=xn+ey/100;
        yn=yn+ex/100;
        e1=xn-subs(eqx);
        e2=yn-subs(eqy);
    end
    
    end

I obtain a good aproximation when i compared to results obtained with commercial softwares for soil mechanics.

I know that my code do not work very well but the aproximation is good.

I can't used fsolve because i don't have that toolbox.

Thanks to everyone.

[John D'Errico]

Note that recursion will find a solution, if it converges. This is called fixed point iteration. If it convergeges to a solution then good, but only if you have found the correct solution.

Why do I say this? Because fixed point iteration need not converge to a solution in some circumstances, and that may be the solution you need to find. That is, fixed point will diverge away from entirely valid solutions, if the derivatives of the right hand side are large near the solution.

The point is, fixed point is a pot luck thing. It may work. It may not. And there is no indication that this will be a problem, because it will simply avoid some valid solutions. As I say, pot luck. You get what you will get.

Answer 1

[John D'Errico]

Assuming all other parameters are known...

Suppose you were to clear all the fractions? So cross multiply by everything in the demoninators? When you got done, you would have two polynomials in the TWO unknown variables. I'm not telling you to do this. But it will teach you something, if you follow along with what I say, in this little thought experiment.

Now collect all of the terms together, as coefficients of one of those polynomials in x_n. The coefficients will be functions of y, but ignore that for the moment. This plynomial appears to be at least a 5th degree polynomial. So there will be an x^5 in that mess, possibly higher order terms too. I'm lazy, and I won't look too hard. The point is, you have what is equivalent to a 5th degree polynomial. If you could find a solution to it, it would look like a really messy version of the quadratic formula, but way more messy. Lots of radicals, etc. And all of that a function of y. But IF you could do this, it would effectively give you 5 roots, each of which was a function of y_n.

Of course, all of that is a mere thought experiment, because you cannot in fact solve for the roots of that polynomial. Abel-Ruffini proved that an infinite number of years ago. (Before I was born, so as far as I am concerned, infinitely long ago.) 5th degree of higher, there are no solutions available in general to the problem. (Some rare polynomials will have solutions, but pretty much never the ones you will have in practice. And not this case for sure.)

Anyway, IF you could solve that, each root will be a function of y. Now plug each root from that mess back into the other equation. It too will be at least quintic in degree, and since roots for x will also be a function of y, the result will now be much higher in degree. As a guess, this might be as much as a 25th degree polynomial for y. It would ONLY involve y now. And unfortunately, there is no way to compute that polynomial, because we are not able to sove the first problem. So we cannot actually perform the mathematics to compute the coefficients of this envisioned polynomial in y.

If you could factor that final polynomial in y, it would have something like 25 roots, many of them complex. It would have a FINITE number of roots. But there would be a lot of them.

Again, you cannot do all of that. But there are some things you CAN do. You can use tools like fimplicit to plot solutions. You CAN use tools like fsolve to find solutions. But fsolve can find only one solutino at a time, and it will depend on where yoyu start the solver, which solution it finds.

Remember though, if ANY of the unknown constants in your probllem are actually unknown? And you want to solve this as a function of the unknown constats too? Now EVERYTHING is now impossible. There is essentially nothing you can do.

I'll try an example, because your problem is way to messy to do here, and I lack any of the unknown parameters anyway.

For example, consider the problem

x = (1 + 1/x) + (1 - 1/y^2) + 3

y = (1 + y/x) + (1 - 1/y) + 2

It is sort of similar to yours, but way simpler. We can first solve one of the equations for x, just to see what will happen, and how things explode.

syms x y
xsol = solve(-x + (1 + 1/x) + (1 - 1/y^2) + 3,'returnconditions',true)

xsol = struct with fields:
             x: [2×1 sym]
    parameters: [1×0 sym]
    conditions: [2×1 sym]

xsol.x

ans = 

xsol.conditions

ans = 

So as long as x and y are not zero, then we have a solution for x, as a function of y. It is pretty messy. But now, we could substitute that solution back into the other equation.

yeq = subs(-y + (1 + y/x) + (1 - 1/y) + 2,x,xsol.x)

yeq = 

The result being two rather messy looking equations for y. Do they have a solution? Even here, if we cleared out the demoninators, cleared the square roots, etc., the result would be at least a 5th degree polynomial in y.

vpasolve(yeq(1))

ans = 

0.62959488914733937908590302617341

A quick try with vpasolve found a solution, but that does not mean those are all of the solutions.

fplot(yeq(1)),yline(0,'r')

In fact, we can see a plot found two solutions for y. Now we could use guess from the plots to find the actual roots, and then go back and solve for x.

All of this was pretty easy, at least relatively so, but here I was able to solve one of the equations for x. Suppose that was not possible, as in your case?

fun1 = @(x,y) -x+ (1 + 1./x) + (1 - 1./y.^2) + 3;
fun2 = @(x,y) -y + (1 + y./x) + (1 - 1./y) + 2;

fimplicit(fun1,'r',[-2,2 -2,2])
hold on
fimplicit(fun2,'b',[-2,2 -2,2])
xlabel x
ylabel y

Wherever the red and blue curves cross, is a solution. Now use those points found graphically, as starting values for fsolve. For example, a solution appears to live near (-.5,.5).

format long g
funs = @(xy) [fun1(xy(1),xy(2));fun2(xy(1),xy(2))];
[xysol,fval,exitflag] = fsolve(funs,[-.5,5])

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

xysol = 1×2

-0.35329172607673 0.629594889147338

fval = 2×1

1.0e+00 * -1.93178806284777e-13 -1.00364161426114e-13

exitflag = 
     1

Now we could try again to find a root in the first quadrant.

[xysol,fval,exitflag] = fsolve(funs,[0.1,0.2])

No solution found.

fsolve stopped because the problem appears regular as measured by the gradient,
but the vector of function values is not near zero as measured by the
value of the function tolerance.

xysol = 1×2

0.999999234604208 0.431789974350947

fval = 2×1

-0.363580058462071 1.68405957363718

exitflag = 
    -2

This has failed. I might need better starting values. Such is life.

Answer 2

[Walter Roberson]

Construct a function that takes a vector of two values, and returns a vector of length 2. The first element of the vector should be the right hand side of the first equation minus the left hand side, and the second element should be the right hand side of the second equation minus the left hand side. Now you can use fsolve()