# how to handle a “Stiff” algebraic equation numerically?

I have a question of great practical importance for me, but I would like to ask it on a bit more of a theoretical mode, because I feel I lack the basic knowledge on it. I would also like to mention that I am a student in physics, not in mathematics, so I apologize in advance if I lag behind on some notions

I want to numerically solve an algebraic equation (I was using Matlab to try and solve it). At some points in the function that I have to define, I stumble upon a term that looks somewhat like :

$$e^{ (A-e^x)B }$$

where "B >> 1" (let's say B=10^6). Let's assume here that A is 1 and x is the unknown parameter I want to solve for.

Here is the problem that I encounter : there is very very few flexibility on the acceptable values of x for matlab.

Let's say we have an initial guess corresponding to a value of 1.1 for exp(x). Then $$(1-e^x)B = 10^5$$, and the term evaluates as $$e^{10^5}$$ which is way beyond what matlab can handle ($$realmax = 1.8e308$$). Matlab understands the term as an infinite value, which, as you can guess, leads to many problems for further solving. Actually, a homebrew implementation of Halley's method shows me that the value of this term switches directly from 0 to inf as (A-exp(x))*B switches sign during iterations ! This is an indication that $$|(A-e^x)B|$$ is simply too large unless x is very close to ln(A)

In fact, A-exp(x) has to be numerically very small for matlab to handle the thing. Which means that the search interval should be extremely narrow and x always very close to ln(A). The solution in general has no reason to be ln(A) (there are many other terms in the eqution apart from this problematic one). Any algorithm that fidgets around a little bit too much will end up having troubles of inf values.

My first question is : how do we describe this class kind of sensitive numerical problems with algebraic equations ? Surely, this is nothing new, and many people more clever than I am have worked on this kind of issue before, so I would like to look into it with the right keywords. The closest idea that I could think of, would be the "stiffness" of ODEs, but here I am working with a purely algebraic system.

Second : what kind of workaround could I try? Of course, the naive one would be to just increase the realmax value of matlab, which should help get my way around these annoying inf values. But I don't think it's possible, and probably not a good idea to do it anyway with numbers in the order of exp(10^6~7). And even if I could do that, I suppose anyway that having a term suddenly switching from 0 to exp(10^6) is a surprise that would be very poorly appreciated by most solvers...

Is there some famous algorithms that I could try to implement ? Thanks in advance !

N.B:

I have tried several algorithms, with different reformulations of my equations, but in the end, it seems it always comes down to the same problem in different forms. If I try to reformulate with log functions for example, the thing in the log turns out to be very sensitive too, and become very easily negative as the algorithm tries different values of x, causing similar problems.

Otherwise, I thought of trying to set bounds for x so that it doesn't go insane. But this turns out to have the same kind problems with a lot of log( -...), very similar to the actual solving. Which makes sense because finding such a precise search interval for x would pretty much mean solving the equation for me :D

EDIT: If it may help, my complete equation may be written under the following form:

$$0=C-\beta V_{sh}(x) - \gamma .exp(\frac{V_{sh}(x)-x}{\delta})$$

where

$$V_{sh}(x)=(A-e^{\alpha x})B$$

with " $$\frac{B}{\delta} >> 1$$ ", and the problematic term is $$exp(\frac{V_{sh}(x)}{\delta})$$

When you have to solve nonlinear equations, the basic idea is to make them as close to linearity as possible in order to solve it in a minimum numer of iterations (remember that Newton method assumes that locally the function is linear).

Consider the case where you want to solve $$e^{ (A-e^{f(x)})B }=k \qquad \text{with} \qquad k >0$$Start taking logarithms to get $$(A-e^{f(x)} )B=\log(k)\implies f(x)=\log \left(A-\frac{\log (k)}{B}\right)$$ Don't you think that this could be much easier ?

• Unfortunately, i do not think so. I am aware that the equation you propose would not require numerical calculations. As mentionned in my nota bene, I already tried such log-based linearization strategies. However, this term is only a part of a more complex equation. The term that will be enclosed in the ln( ) in the end will depend on x in a non-linear way. It will be very likely to become negative as the algorithm fidgets on x, in a similar manner as the exponential was likely to become infinite. This reformulation ends up producing just as many errors – Barbaud Julien Jan 17 at 6:50
• @BarbaudJulien. Could you provide an example ? If you prefer, my e-mail address is in my profile. By the way, where are you in France ? I am in Pau. – Claude Leibovici Jan 17 at 6:54
• I have included the complete equation in my question as an edit rather than a mail, because I think it can be beneficial if someone else wants to give it a shot. And I am from France, but right now living in Shanghai ^^ – Barbaud Julien Jan 17 at 8:03
• @BarbaudJulien. Could you provide a set of values for $(\alpha,\beta,\gamma,\delta,A,B)$ of your choice ? I would really like to play with this equation. Cheers :-) – Claude Leibovici Jan 17 at 9:47