Re: vpasolve to find all real roots for multiple nonlinear equations and multiple variables

[Oliver Heimlich]

On 01.06.2016 08:38, sam kjm wrote:
> Thanks Oliver for the guidance! I think I managed  to get some results.
> However I met some problems along the way again.
> 
> 
> I have 2 separate code files:
> 
> 1) My function file ("steady5.m")
> 
> function steadystate = steady5(y, total_prot2)
> 

...

> 
>     ## Unknowns to be solved (prot1, prot2s)
>     prot1 = y(1)
>     prot2s = y(2)
>     prot2 = total_prot2 - prot2s;

...

> end
> 
> 

> 3) I got the results i think. But fsolve displays a lot of iterations before
> I get the final solution through disp(y). Unless I actually did not get the
> solutions? 
> 
> Plus what if I want to get >2 solutions? I have some that should have 3
> solutions but I didn't manage to replicate that.  Is it determined by the
> initial values infsup([], [], []) too?
> 

I have replayed your example. A lot of output comes from the missing
semicolons in steady5. The fsolve function bisects the search range and
calls steady5 a lot of times (brute-force method). That's why you see so
much noise. The final solution should be correct, but it is only an
enclosure of many solutions.

total_protein2 = 0
   [0.23437, 2.002]
       [0, 0.88868]
total_protein2 = 0.50000
   [0.23437, 0.30274]
    [0.26367, 1.0547]

You can get a detailed look on many solutions like this in call_steady5.m:

clear;clc;clf;
total_prot2 = 0;

while total_prot2 < 1
    disp(total_prot2);
    [y, solutions, inner] = fsolve(@(y) steady5(y, total_prot2), infsup
([0;0], [5;5]));

    ## Outer approximation of all solutions
    disp(y);

    clf
    hold on
    ##  Outer approximation of solution set (green)
    plot (solutions(1, :), solutions(2, :))
    ##  Inner approximation of solution set (red)
    plot (solutions(1, inner), solutions(2, inner), [1 0 0])
    xlabel ("prot1")
    ylabel ("prot2s")

    input ("Continue? Press Enter.");

    total_prot2 = total_prot2 + 0.5;
end

As you can see in the plot, there are only small ranges where solutions
might occur. Any white area in the plot is guaranteed to have no
solution. However, the algorithm did not find an inner approximation of
the solution set. That means, you have no guarantee that the green areas
contain any solutions at all. Please keep that in mind!

For example, when you increase the precision & runtime with
optimset("TolX", 0.001, "MaxIter", 25), you will see that the solution
set is actually even smaller. It looks like there were two solutions,
but I might be wrong.

I could reduce the overestimation (and increase the speed) of the
algorithm further by using the following rate expressions in steady5.m.
As a rule of thumb you should try to have each interval variable only
once per expression. Otherwise you'd introduce the so called “dependency
error”, a common problem in interval arithmetic.

    ## setting the rate expressions (degradation rate, production rate etc)
    v0 = k0;
    v1 = k1 / (j1 / prot2 + 1);
    vm1 = km1 / (jm1 / prot2s + 1);
    v2 = (k2 * prot1) / (j2 / prot2s + 1);
    vm3 = (km3 * prot1) / (jm3 / prot2s + 1);

I hope this helps.

Oliver


P.S. In your computations you assume that your model and fixed
parameters are exact. You may also play with inexact parameters and see
what happens. For example, you can use k2 = infsup ("0.4?") to add an
uncertainty of +/- 0.5 decimal ULP to that value.