Thread Subject:

Solving an equation with no close analytical solution

Hi to all and thanks in advance for your help.
I have to solve an equation that writes as
a=exp(x*c)*(exp((b-x)*c)-1)/(b-x)
The unknown is x while a,b and c are known
I have 5000 values for a while c and b are the same through each equation, this means that I need to find 5000 x.
In maple I can find a numerical solution with fsolve, but I have to do manually for each of the 5000 values (which takes a awful amount of time)
I would like to know if there is a way to find a numerical solution to this equation in matlab and to this in a loop that allows to use the array of 5000 values of a to create an array for the 5000 values of x.

I do know if that's clear, i have been struggling for this porblem for the latest three days.
So any help would be really greatly appreciated!!!
thanks again
a

"deva " <andadodk@yahoo.fr> wrote in message <gifvtb$6tp$1@fred.mathworks.com>...
> Hi to all and thanks in advance for your help.
> I have to solve an equation that writes as
> a=exp(x*c)*(exp((b-x)*c)-1)/(b-x)
> The unknown is x while a,b and c are known
> I have 5000 values for a while c and b are the same through each equation, this means that I need to find 5000 x.
> In maple I can find a numerical solution with fsolve, but I have to do manually for each of the 5000 values (which takes a awful amount of time)
> I would like to know if there is a way to find a numerical solution to this equation in matlab and to this in a loop that allows to use the array of 5000 values of a to create an array for the 5000 values of x.
>
> I do know if that's clear, i have been struggling for this porblem for the latest three days.
> So any help would be really greatly appreciated!!!
> thanks again
> a

The simple answer is a loop using fzero, not fsolve.
Fzero will be more robust for problems of one variable,
if you are able to bracket the root.

The more complex answer is to use arrayfun, which
will solve the problem using an implicit loop. Still use
fzero though.

My next choice as an answer is to solve it as a batched
problem in fsolve. Use the largescale method. Specify
the pattern of the jacobian matrix through optimset.
I discuss this style of solution in my optimization tips
and tricks.

http://www.mathworks.com/matlabcentral/fileexchange/8553

You may need to verify the solutions for convergence,
since some subproblems may possibly diverge. Go back
and apply fzero to those subproblems only.

HTH,
John

"deva " <andadodk@yahoo.fr> wrote in message <gifvtb$6tp$1@fred.mathworks.com>...
> Hi to all and thanks in advance for your help.
> I have to solve an equation that writes as
> a=exp(x*c)*(exp((b-x)*c)-1)/(b-x)
> The unknown is x while a,b and c are known
> I have 5000 values for a while c and b are the same through each equation, this means that I need to find 5000 x.
> In maple I can find a numerical solution with fsolve, but I have to do manually for each of the 5000 values (which takes a awful amount of time)
> I would like to know if there is a way to find a numerical solution to this equation in matlab and to this in a loop that allows to use the array of 5000 values of a to create an array for the 5000 values of x.
>
> I do know if that's clear, i have been struggling for this porblem for the latest three days.
> So any help would be really greatly appreciated!!!
> thanks again
> a

  My symbolic toolbox finds a solution in terms of the LambertW function:

 f = 'exp(c*x)*(exp((b-x)*c)-1)/(b-x)-a'
 x = solve('exp(c*x)*(exp((b-x)*c)-1)/(b-x)-a','x')
 x =
 -(lambertw(-exp((-exp(c*b)+a*b)*c/a)/a*c)+1/a*c*exp(c*b)-c*b)/c

I notice further that the file exchange has a LambertW function written by Pascal Getreuer at

 http://www.mathworks.com/matlabcentral/fileexchange/6909

It even states that it handles vector inputs, so you could write a single vectorized expression to solve your problem for the 5000 values of 'a' all in one line. I will confess I have never used this function so you can tell us how well it works for you.

Roger Stafford

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gihaei$pim$1@fred.mathworks.com>...
> "deva " <andadodk@yahoo.fr> wrote in message <gifvtb$6tp$1@fred.mathworks.com>...
> > Hi to all and thanks in advance for your help.
> > I have to solve an equation that writes as
> > a=exp(x*c)*(exp((b-x)*c)-1)/(b-x)
> > The unknown is x while a,b and c are known
> > I have 5000 values for a while c and b are the same through each equation, this means that I need to find 5000 x.
> > In maple I can find a numerical solution with fsolve, but I have to do manually for each of the 5000 values (which takes a awful amount of time)
> > I would like to know if there is a way to find a numerical solution to this equation in matlab and to this in a loop that allows to use the array of 5000 values of a to create an array for the 5000 values of x.
> >
> > I do know if that's clear, i have been struggling for this porblem for the latest three days.
> > So any help would be really greatly appreciated!!!
> > thanks again
> > a
>
> My symbolic toolbox finds a solution in terms of the LambertW function:
>
> f = 'exp(c*x)*(exp((b-x)*c)-1)/(b-x)-a'
> x = solve('exp(c*x)*(exp((b-x)*c)-1)/(b-x)-a','x')
> x =
> -(lambertw(-exp((-exp(c*b)+a*b)*c/a)/a*c)+1/a*c*exp(c*b)-c*b)/c
>
> I notice further that the file exchange has a LambertW function written by Pascal Getreuer at
>
> http://www.mathworks.com/matlabcentral/fileexchange/6909
>
> It even states that it handles vector inputs, so you could write a single vectorized expression to solve your problem for the 5000 values of 'a' all in one line. I will confess I have never used this function so you can tell us how well it works for you.
>
> Roger Stafford
Hi,
thanks for your message.
I admit that I found the solution with the lambertw when I first tried to solve he equation, however the way it works (at least this seems to me) is that to find solution x comes out equal to 0.0017 that is exactly the value of b.
So at some point to get this solution the solver assumes b=x, that is not what I want.
When I try to solve it numerically by substitution to obtain the value of x I get 0.005835 by putting a=5007, b=0.0017 and c=540 THe value of 0.005835 is what the solver in maple (fsolve) gives to me.
But I cannot do the same thing in matlab, even not with the fzero solver....
I don't know why
thanks a lot again

"John D'Errico" <woodchips@rochester.rr.com> wrote in message <gigbvi$pfg$1@fred.mathworks.com>...
> "deva " <andadodk@yahoo.fr> wrote in message <gifvtb$6tp$1@fred.mathworks.com>...
> > Hi to all and thanks in advance for your help.
> > I have to solve an equation that writes as
> > a=exp(x*c)*(exp((b-x)*c)-1)/(b-x)
> > The unknown is x while a,b and c are known
> > I have 5000 values for a while c and b are the same through each equation, this means that I need to find 5000 x.
> > In maple I can find a numerical solution with fsolve, but I have to do manually for each of the 5000 values (which takes a awful amount of time)
> > I would like to know if there is a way to find a numerical solution to this equation in matlab and to this in a loop that allows to use the array of 5000 values of a to create an array for the 5000 values of x.
> >
> > I do know if that's clear, i have been struggling for this porblem for the latest three days.
> > So any help would be really greatly appreciated!!!
> > thanks again
> > a
>
> The simple answer is a loop using fzero, not fsolve.
> Fzero will be more robust for problems of one variable,
> if you are able to bracket the root.
>
> The more complex answer is to use arrayfun, which
> will solve the problem using an implicit loop. Still use
> fzero though.
>
> My next choice as an answer is to solve it as a batched
> problem in fsolve. Use the largescale method. Specify
> the pattern of the jacobian matrix through optimset.
> I discuss this style of solution in my optimization tips
> and tricks.
>
> http://www.mathworks.com/matlabcentral/fileexchange/8553
>
> You may need to verify the solutions for convergence,
> since some subproblems may possibly diverge. Go back
> and apply fzero to those subproblems only.
>
> HTH,
> John

Hi John,
thanks a lot for your message.
I have tried and I get the message
"??? Error using ==> fzero
FZERO cannot continue because user supplied function_handle ==> @(x) myfun(x,a,b,c)
failed with the error below.

Error: File: myfun.m Line: 5 Column: 1
Function definitions are not permitted at the prompt or in scripts. "

What I have done is
function l=myfun(x,a,b,c)
l=(c-(exp(x*a).*(exp((b-x)*a)-1))/(b-x));

a=540;
b=(1+0.02)^(1/12)-1;
c=5007;

 fzero(@(x) myfun(x,a,b,c),0.1);

thanks again

"deva " <andadodk@yahoo.fr> wrote in message <giihon$t6b$1@fred.mathworks.com>...

> Hi John,
> thanks a lot for your message.
> I have tried and I get the message
> "??? Error using ==> fzero
> FZERO cannot continue because user supplied function_handle ==> @(x) myfun(x,a,b,c)
> failed with the error below.
>
> Error: File: myfun.m Line: 5 Column: 1
> Function definitions are not permitted at the prompt or in scripts. "
>
> What I have done is
> function l=myfun(x,a,b,c)
> l=(c-(exp(x*a).*(exp((b-x)*a)-1))/(b-x));
>
> a=540;
> b=(1+0.02)^(1/12)-1;
> c=5007;
>
> fzero(@(x) myfun(x,a,b,c),0.1);
>
> thanks again

Read the error message. you cannot just type in
a function definition as above at the command
line. And a script file is equivalent to typing in
the lines contained at the command line.

You can save myfun as a SEPARATE m-file on
your search path. Or you can write myfun as an
anonymous function.

a=540;
b=(1+0.02)^(1/12)-1;
c=5007;
myfun = @(x) (c-(exp(x*a).*(exp((b-x)*a)-1))/(b-x));
 
x = fzero(myfun,0.1)
x =
    0.0058388

Better anyway is to use Roger's suggestion,
writing the solution using the Lambert W.

John

"deva " <andadodk@yahoo.fr> wrote in message <giihj1$ije$1@fred.mathworks.com>...
> .......
> I admit that I found the solution with the lambertw when I first tried to solve he equation, however the way it works (at least this seems to me) is that to find solution x comes out equal to 0.0017 that is exactly the value of b.
> So at some point to get this solution the solver assumes b=x, that is not what I want.
> When I try to solve it numerically by substitution to obtain the value of x I get 0.005835 by putting a=5007, b=0.0017 and c=540 THe value of 0.005835 is what the solver in maple (fsolve) gives to me.
> But I cannot do the same thing in matlab, even not with the fzero solver....
> ........

  Deva, note that the Lambert W function extended to the complex plane has more than one branch. As a consequence, for some real values there are two real solutions. Though I haven't done any computing with your data it sounds from your description as though you might be running into that problem.

  Don't forget that you can always test any answer for accuracy by putting it back into your original equation. You didn't state whether the answer you considered incorrect did or did not satisfy this equation.

  You would also have the same difficulty using matlab's 'fzero'. Fzero's final answer would in that case depend on what initial estimate you furnished it with.

  Pascal Getreuer's matlab function has the option of entering a "branch" argument input which is called 'b', (not to be confused with your 'b'.) I don't know the details but you might explore that possibility as a way of leading to the answer that you consider correct.

  You can consult the Wikipedia website for some more information about the Lambert W function at:

 http://en.wikipedia.org/wiki/Lambert_W_function

Roger Stafford

"deva " <andadodk@yahoo.fr> wrote in message <giihj1$ije$1@fred.mathworks.com>...
> I admit that I found the solution with the lambertw when I first tried to solve he equation, however the way it works (at least this seems to me) is that to find solution x comes out equal to 0.0017 that is exactly the value of b.
> So at some point to get this solution the solver assumes b=x, that is not what I want.
> .........

  I looked a little more closely at your problem, Deva. By multiplying both sides of your equation by b-x you get

 a*(b-x) = exp(x*c)*(exp((b-x)*c)-1) = exp(b*c) - exp(x*c)

and with this form of the equation, it is obvious that x = b is a solution. In this loose sense, this is a solution to your original equation as well as the other value you mentioned. As it turns out, the x = b value is on what is regarded as the "principal" branch of the Lambert W function, so that is what it yielded as a solution. As I mentioned earlier, to get around this you will have to use Pascal Getreuer's second input argument which is called 'b' to force 'lambertw' onto another branch. I do not know what value of 'b' this will require, but I suspect it is b = 1 and if not that then b = -1.

Roger Stafford

"Roger Stafford" <ellieandrogerxyzzy@mindspring.com.invalid> wrote in message <gikasr$a4q$1@fred.mathworks.com>...
> ........
> As I mentioned earlier, to get around this you will have to use Pascal Getreuer's second input argument which is called 'b' to force 'lambertw' onto another branch. I do not know what value of 'b' this will require, but I suspect it is b = 1 and if not that then b = -1.
> ........

  I just tried out 'lambertw'. For your stuff you need to set the "branch" input equal to -1 in order to stay off the principal branch.

  Its output may have a very tiny imaginary part due to round off error and you can take its real part to get rid of that:

  w = real(lambertw(-1,z));

where z is all the stuff that should go inside the Lambert W function.

  However, you should always check that your results satisfy the original equation. It is possible to have 'a' set so large that for the given 'b' and 'c' values there are no real-valued solutions except the one at x = b.

  As I mentioned earlier, you can use a vector form of z for your 5000 points and avoid for-loops altogether.

Roger Stafford