Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
fsolve to nonlinear equation

Subject: fsolve to nonlinear equation

From: edward kabanyas

Date: 14 Apr, 2009 10:04:01

Message: 1 of 7

Dear All,

I would like to find the parameters xx(1) and xx(2) of the following equation:

function f = myfun(xx,eta,Dmax, Dmin)
f = ((gamma(xx(1)+4)*(gammainc(xx(2)*Dmax,xx(1)+4,'lower')-(1-gammainc(xx(2)*Dmin,xx(1)+4, 'upper'))))^2./((gamma(xx(1)+5)*(gammainc(xx(2)*Dmax,xx(1)+5,'lower')-(1-gammainc(xx(2)*Dmin,xx(1)+5, 'upper'))))*(gamma(xx(1)+3)*(gammainc(xx(2)*Dmax,xx(1)+3,'lower')-(1-gammainc(xx(2)*Dmin,xx(1)+3, 'upper'))))))-eta;

eta, Dmax, Dmin are constant and given and initial values of xx are also given. I tried several matlab optimization such as fzero and fsolve. but it does not work. I think fzero does not work because the equation contain 2 variables, but why fsolve does not work too? The following errors occur when I run my code:

 xval = fsolve(@(xx)myfun(xx,eta,Dmax, Dmin),xx);

xval =

    7.7944 14.6764


eta =

    0.8466

Warning: Default trust-region dogleg method of FSOLVE cannot
 handle non-square systems; switching to Gauss-Newton method.
> In fsolve at 232
  In msearch at 122
Optimization terminated: directional derivative along
 search direction less than TolFun and infinity-norm of
 gradient less than 10*(TolFun+TolX).

Any help is really appreciated.

Regards;
Edward
 

Subject: fsolve to nonlinear equation

From: John D'Errico

Date: 14 Apr, 2009 10:23:02

Message: 2 of 7

"edward kabanyas" <djuky_hmi@yahoo.com> wrote in message <gs1n2h$2i$1@fred.mathworks.com>...
> Dear All,
>
> I would like to find the parameters xx(1) and xx(2) of the following equation:
>
> function f = myfun(xx,eta,Dmax, Dmin)
> f = ((gamma(xx(1)+4)*(gammainc(xx(2)*Dmax,xx(1)+4,'lower')-(1-gammainc(xx(2)*Dmin,xx(1)+4, 'upper'))))^2./((gamma(xx(1)+5)*(gammainc(xx(2)*Dmax,xx(1)+5,'lower')-(1-gammainc(xx(2)*Dmin,xx(1)+5, 'upper'))))*(gamma(xx(1)+3)*(gammainc(xx(2)*Dmax,xx(1)+3,'lower')-(1-gammainc(xx(2)*Dmin,xx(1)+3, 'upper'))))))-eta;
>
> eta, Dmax, Dmin are constant and given and initial values of xx are also given. I tried several matlab optimization such as fzero and fsolve. but it does not work. I think fzero does not work because the equation contain 2 variables, but why fsolve does not work too? The following errors occur when I run my code:
>
> xval = fsolve(@(xx)myfun(xx,eta,Dmax, Dmin),xx);
>
> xval =
>
> 7.7944 14.6764
>
>
> eta =
>
> 0.8466
>
> Warning: Default trust-region dogleg method of FSOLVE cannot
> handle non-square systems; switching to Gauss-Newton method.
> > In fsolve at 232
> In msearch at 122

No. This was a WARNING. NOT an error. It
indicated that you allowed the default method
to be used for a problem that could not employ
that method, so it switched methods.

This is all it said. NOT an error.


> Optimization terminated: directional derivative along
> search direction less than TolFun and infinity-norm of
> gradient less than 10*(TolFun+TolX).

Again, NOT an ERROR.

Look in the dictionary. "terminated" will appear as
a synonym for "done". It has finished its work on
this problem. Nothing other than that. Then it
told you why it thought that it had reasonable
grounds for deciding it had finished.

John

Subject: fsolve to nonlinear equation

From: edward kabanyas

Date: 14 Apr, 2009 10:34:02

Message: 3 of 7

Dear John;

thanks for your message. I think it reflects error because no change from initial value to the output one.

Again, thanks.
Edward

"John D'Errico" <woodchips@rochester.rr.com> wrote in message <gs1o66$9t2$1@fred.mathworks.com>...
> "edward kabanyas" <djuky_hmi@yahoo.com> wrote in message <gs1n2h$2i$1@fred.mathworks.com>...
> > Dear All,
> >
> > I would like to find the parameters xx(1) and xx(2) of the following equation:
> >
> > function f = myfun(xx,eta,Dmax, Dmin)
> > f = ((gamma(xx(1)+4)*(gammainc(xx(2)*Dmax,xx(1)+4,'lower')-(1-gammainc(xx(2)*Dmin,xx(1)+4, 'upper'))))^2./((gamma(xx(1)+5)*(gammainc(xx(2)*Dmax,xx(1)+5,'lower')-(1-gammainc(xx(2)*Dmin,xx(1)+5, 'upper'))))*(gamma(xx(1)+3)*(gammainc(xx(2)*Dmax,xx(1)+3,'lower')-(1-gammainc(xx(2)*Dmin,xx(1)+3, 'upper'))))))-eta;
> >
> > eta, Dmax, Dmin are constant and given and initial values of xx are also given. I tried several matlab optimization such as fzero and fsolve. but it does not work. I think fzero does not work because the equation contain 2 variables, but why fsolve does not work too? The following errors occur when I run my code:
> >
> > xval = fsolve(@(xx)myfun(xx,eta,Dmax, Dmin),xx);
> >
> > xval =
> >
> > 7.7944 14.6764
> >
> >
> > eta =
> >
> > 0.8466
> >
> > Warning: Default trust-region dogleg method of FSOLVE cannot
> > handle non-square systems; switching to Gauss-Newton method.
> > > In fsolve at 232
> > In msearch at 122
>
> No. This was a WARNING. NOT an error. It
> indicated that you allowed the default method
> to be used for a problem that could not employ
> that method, so it switched methods.
>
> This is all it said. NOT an error.
>
>
> > Optimization terminated: directional derivative along
> > search direction less than TolFun and infinity-norm of
> > gradient less than 10*(TolFun+TolX).
>
> Again, NOT an ERROR.
>
> Look in the dictionary. "terminated" will appear as
> a synonym for "done". It has finished its work on
> this problem. Nothing other than that. Then it
> told you why it thought that it had reasonable
> grounds for deciding it had finished.
>
> John

Subject: fsolve to nonlinear equation

From: John D'Errico

Date: 14 Apr, 2009 11:13:01

Message: 4 of 7

"edward kabanyas" <djuky_hmi@yahoo.com> wrote in message <gs1oqq$kv3$1@fred.mathworks.com>...
> Dear John;
>
> thanks for your message. I think it reflects error because no change from initial value to the output one.
>

Do you think you might have said this the first time?

Again, this is not an error message. It is a simple
statement of operation. The algorithm has attempted
to iterate to a solution. If you are not happy with
where it has gone (or apparently did not go) that is
another matter.

Lack of progress made is a probable indication of
user error. Often the user has posed an objective
with a bug in it, or passed the function or parameters
improperly. Another VERY common reason for failure
of this sort is a terrible choice of starting values. (I
think this to be one actual reason fsolve is failing to
go anywhere.)

I cannot know what your problem is, because I
cannot test your code. You have told us the value
of eta here, but not Dmax or Dmin. I assume the
starting value you chose was the one that was
returned from your comments. How did you choose
the starting values?

The next problem is that it appears from reading
your code that you have a problem with ONE
equation in TWO unknowns. Assuming that a
solution does exist, this is not really a problem
for fsolve. If a solution does exist, then there will
probably be infinitely many solutions. In fact,
I would suggest using a simple contour plot to
display the solution locus. Display only the
contour at zero.

I cannot do this myself, since you did not give
enough information. (Dmin = ? Dmax = ?)

john

Subject: fsolve to nonlinear equation

From: edward kabanyas

Date: 14 Apr, 2009 12:46:01

Message: 5 of 7

Dear John;

Thanks again for your message. The following is the example parameter being used:

Dmax = 8;
Dmin = 0.05;
eta = 0.8392 (calculated from measured data);

Initial values of xx = [2.2190, 6.1253]

The result of xval = fsolve(@(xx)myfun(xx,eta,Dmax, Dmin),xx) is

2.2190 6.1253

So we can see no change between input and output file. Initial value of x are calculated using gamma function (untrancated). But I need the values of x for truncated gamma (incomplete gamma function) given by:

f = ((gamma(xx(1)+4)*(gammainc(xx(2)*Dmax,xx(1)+4,'lower')-(1-gammainc(xx(2)*Dmin,xx(1)+4, 'upper'))))^2./((gamma(xx(1)+5)*(gammainc(xx(2)*Dmax,xx(1)+5,'lower')-(1-gammainc(xx(2)*Dmin,xx(1)+5, 'upper'))))*(gamma(xx(1)+3)*(gammainc(xx(2)*Dmax,xx(1)+3,'lower')-(1-gammainc(xx(2)*Dmin,xx(1)+3, 'upper'))))))-eta;

For gamma (untrancated) we can solve f numerically, so we can get xx(1) and xx(2) easily, but for incomplete gamma function, I should use iterative procedure to find xx(1) and xx(2) by using the values from untrancated gamma as initial.

Actually the difference between truncated and untrancated is about max 50%
and parameter of truncated < untrancated parameters.

Hope you could get the point. Again, thanks very much.

Regards
Edward

Subject: fsolve to nonlinear equation

From: John D'Errico

Date: 14 Apr, 2009 14:34:23

Message: 6 of 7

"edward kabanyas" <djuky_hmi@yahoo.com> wrote in message <gs20i9$ea6$1@fred.mathworks.com>...
> Dear John;
>
> Thanks again for your message. The following is the example parameter being used:
>
> Dmax = 8;
> Dmin = 0.05;
> eta = 0.8392 (calculated from measured data);
>
> Initial values of xx = [2.2190, 6.1253]
>
> The result of xval = fsolve(@(xx)myfun(xx,eta,Dmax, Dmin),xx) is
>
> 2.2190 6.1253
>
> So we can see no change between input and output file. Initial value of x are calculated using gamma function (untrancated). But I need the values of x for truncated gamma (incomplete gamma function) given by:
>
> f = ((gamma(xx(1)+4)*(gammainc(xx(2)*Dmax,xx(1)+4,'lower')-(1-gammainc(xx(2)*Dmin,xx(1)+4, 'upper'))))^2./((gamma(xx(1)+5)*(gammainc(xx(2)*Dmax,xx(1)+5,'lower')-(1-gammainc(xx(2)*Dmin,xx(1)+5, 'upper'))))*(gamma(xx(1)+3)*(gammainc(xx(2)*Dmax,xx(1)+3,'lower')-(1-gammainc(xx(2)*Dmin,xx(1)+3, 'upper'))))))-eta;
>
> For gamma (untrancated) we can solve f numerically, so we can get xx(1) and xx(2) easily, but for incomplete gamma function, I should use iterative procedure to find xx(1) and xx(2) by using the values from untrancated gamma as initial.
>
> Actually the difference between truncated and untrancated is about max 50%
> and parameter of truncated < untrancated parameters.
>
> Hope you could get the point. Again, thanks very much.

You have a SINGLE equation.

THERE ARE INFINITELY MANY SOLUTIONS.

Read my last response. AGAIN. This is not a problem
for fsolve. Since you have now provided enough
information to look at the function itself, lets try
it out. See that I will NEVER call fsolve, certainly for
a long time. Look at what is happening.

>> Dmax = 8;
>> Dmin = 0.05;
>> eta = 0.8392;
>> xx0 = [2.2190, 6.1253];

I've defined the function as an m-file,
called myfun.m. The VERY first thing I will
do is to evaluate the function at xx0. Do
not do ANYTHING else before you do this.
Look at what is happening. THINK about
what you see.

>> myfun(xx0,eta,Dmax,Dmin)
ans =
   9.4193e-06

So the very first thing I see is that this is a
SCALAR function of TWO variables. As
important, the function value is already a
small number, quite near zero in fact.

Try different points near the starting value.

>> myfun([2 6],eta,Dmax,Dmin)
ans =
   -0.0058548
>> myfun([2.3 6.4],eta,Dmax,Dmin)
ans =
    0.0020767

So, while the function is not constant, the
starting value is already a root to within the
default tolerance that fsolve will use. I will not
yet even think about calling fsolve though.

UNDERSTAND what you are doing. THINK about
your function, and what it looks like. How can
we do this most easily? By plotting it. It is just
a function of the form z(x,y).

>> [x1,x2] = meshgrid(2:.01:2.4,6:.01:6.2);
>> fv = zeros(size(x1));
>> for i = 1:numel(x1)
>> fv(i) = myfun([x1(i),x2(i)],eta,Dmin,Dmax);
>> end

Plot the function. Understand what it looks like,
before you blindly throw a root finder at it. This
is especially valuable when you are having
problems if you DID already throw a root finder
at it and did not have success.

PLOT THE FUNCTION. Just do it!

>> surf(x1,x2,fv)

No surprise. The function is pretty much a planar
surface. The "root" is not a single root at all, but
an infinite number of roots, all of which will lie
along a single (possibly slightly curved) line in the
parameter space.

How can we look at this line, this locus of solutions
to the single equation in two unknowns? Do it
exactly as I told you to do so in my response. Use
contour. Since you have already subtracted off the
value of eta in myfun, just plot the zero contour.

>> contour(x1,x2,fv,[0 0])
>> xlabel x1
>> ylabel y1
>> grid on

See that the solution locus is a vertical line, with
x1 essentially fixed at one level. There is still NO
reason to use fsolve yet.

The contour plot told us the solution locus is
essentially a straight line. I'll use my derivest
function to look at the shape of the function
in this neighborhood.

>> format long g
>> derivest(@(xx1) myfun([xx1,xx0(2)],eta,Dmin,Dmax),xx0(1),'style','forward','vect','no')
ans =
     0.0258357008522912

myfun([xx0(1),xx2],eta,Dmin,Dmax),xx0(2),'style','forward','vect','no')
ans =
      5.53545630592774e-06

>> myfun(xx0,eta,Dmin,Dmax)
ans =
      9.41932297959802e-06

So a VERY good LINEAR model for this function is

z = 9.41932297959802e-06 + 0.0258357008522912*(x1 - 2.219) + 5.53545630592774e-06*(x2 - 6.1253)

What is THE root of this function? Do you appreciate
that it has infinitely many roots, and that they fall along
a line?

Pick some round value for x2, say 6.1, and solve for
x1, such that z = 0. Let me see how good the model is.

>> x1 = 2.219 + (- 5.53545630592774e-06*(6.1 - 6.1253) - ...
    9.41932297959802e-06)/0.0258357008522912

x1 =
          2.21864083513785
 
Now, let me try it out in myfun.

>> myfun([x1,6.1],eta,Dmin,Dmax)
ans =
      4.24534185583525e-10

Indeed, this is pretty close to zero.

Have I ever touched fsolve? NO. Look at your
problem. Think about it. Don't just blindly assume
the computer will be smart enough to solve your
problem with no thought of yours invested. This
only works in the movies or on some mindless TV
show.

John

Subject: fsolve to nonlinear equation

From: edward kabanyas

Date: 14 Apr, 2009 14:58:01

Message: 7 of 7

Dear John;

Thanks for your suggestion.

Regards;
Edward
"John D'Errico" <woodchips@rochester.rr.com> wrote in message <gs26tf$ph2$1@fred.mathworks.com>...
> "edward kabanyas" <djuky_hmi@yahoo.com> wrote in message <gs20i9$ea6$1@fred.mathworks.com>...
> > Dear John;
> >
> > Thanks again for your message. The following is the example parameter being used:
> >
> > Dmax = 8;
> > Dmin = 0.05;
> > eta = 0.8392 (calculated from measured data);
> >
> > Initial values of xx = [2.2190, 6.1253]
> >
> > The result of xval = fsolve(@(xx)myfun(xx,eta,Dmax, Dmin),xx) is
> >
> > 2.2190 6.1253
> >
> > So we can see no change between input and output file. Initial value of x are calculated using gamma function (untrancated). But I need the values of x for truncated gamma (incomplete gamma function) given by:
> >
> > f = ((gamma(xx(1)+4)*(gammainc(xx(2)*Dmax,xx(1)+4,'lower')-(1-gammainc(xx(2)*Dmin,xx(1)+4, 'upper'))))^2./((gamma(xx(1)+5)*(gammainc(xx(2)*Dmax,xx(1)+5,'lower')-(1-gammainc(xx(2)*Dmin,xx(1)+5, 'upper'))))*(gamma(xx(1)+3)*(gammainc(xx(2)*Dmax,xx(1)+3,'lower')-(1-gammainc(xx(2)*Dmin,xx(1)+3, 'upper'))))))-eta;
> >
> > For gamma (untrancated) we can solve f numerically, so we can get xx(1) and xx(2) easily, but for incomplete gamma function, I should use iterative procedure to find xx(1) and xx(2) by using the values from untrancated gamma as initial.
> >
> > Actually the difference between truncated and untrancated is about max 50%
> > and parameter of truncated < untrancated parameters.
> >
> > Hope you could get the point. Again, thanks very much.
>
> You have a SINGLE equation.
>
> THERE ARE INFINITELY MANY SOLUTIONS.
>
> Read my last response. AGAIN. This is not a problem
> for fsolve. Since you have now provided enough
> information to look at the function itself, lets try
> it out. See that I will NEVER call fsolve, certainly for
> a long time. Look at what is happening.
>
> >> Dmax = 8;
> >> Dmin = 0.05;
> >> eta = 0.8392;
> >> xx0 = [2.2190, 6.1253];
>
> I've defined the function as an m-file,
> called myfun.m. The VERY first thing I will
> do is to evaluate the function at xx0. Do
> not do ANYTHING else before you do this.
> Look at what is happening. THINK about
> what you see.
>
> >> myfun(xx0,eta,Dmax,Dmin)
> ans =
> 9.4193e-06
>
> So the very first thing I see is that this is a
> SCALAR function of TWO variables. As
> important, the function value is already a
> small number, quite near zero in fact.
>
> Try different points near the starting value.
>
> >> myfun([2 6],eta,Dmax,Dmin)
> ans =
> -0.0058548
> >> myfun([2.3 6.4],eta,Dmax,Dmin)
> ans =
> 0.0020767
>
> So, while the function is not constant, the
> starting value is already a root to within the
> default tolerance that fsolve will use. I will not
> yet even think about calling fsolve though.
>
> UNDERSTAND what you are doing. THINK about
> your function, and what it looks like. How can
> we do this most easily? By plotting it. It is just
> a function of the form z(x,y).
>
> >> [x1,x2] = meshgrid(2:.01:2.4,6:.01:6.2);
> >> fv = zeros(size(x1));
> >> for i = 1:numel(x1)
> >> fv(i) = myfun([x1(i),x2(i)],eta,Dmin,Dmax);
> >> end
>
> Plot the function. Understand what it looks like,
> before you blindly throw a root finder at it. This
> is especially valuable when you are having
> problems if you DID already throw a root finder
> at it and did not have success.
>
> PLOT THE FUNCTION. Just do it!
>
> >> surf(x1,x2,fv)
>
> No surprise. The function is pretty much a planar
> surface. The "root" is not a single root at all, but
> an infinite number of roots, all of which will lie
> along a single (possibly slightly curved) line in the
> parameter space.
>
> How can we look at this line, this locus of solutions
> to the single equation in two unknowns? Do it
> exactly as I told you to do so in my response. Use
> contour. Since you have already subtracted off the
> value of eta in myfun, just plot the zero contour.
>
> >> contour(x1,x2,fv,[0 0])
> >> xlabel x1
> >> ylabel y1
> >> grid on
>
> See that the solution locus is a vertical line, with
> x1 essentially fixed at one level. There is still NO
> reason to use fsolve yet.
>
> The contour plot told us the solution locus is
> essentially a straight line. I'll use my derivest
> function to look at the shape of the function
> in this neighborhood.
>
> >> format long g
> >> derivest(@(xx1) myfun([xx1,xx0(2)],eta,Dmin,Dmax),xx0(1),'style','forward','vect','no')
> ans =
> 0.0258357008522912
>
> myfun([xx0(1),xx2],eta,Dmin,Dmax),xx0(2),'style','forward','vect','no')
> ans =
> 5.53545630592774e-06
>
> >> myfun(xx0,eta,Dmin,Dmax)
> ans =
> 9.41932297959802e-06
>
> So a VERY good LINEAR model for this function is
>
> z = 9.41932297959802e-06 + 0.0258357008522912*(x1 - 2.219) + 5.53545630592774e-06*(x2 - 6.1253)
>
> What is THE root of this function? Do you appreciate
> that it has infinitely many roots, and that they fall along
> a line?
>
> Pick some round value for x2, say 6.1, and solve for
> x1, such that z = 0. Let me see how good the model is.
>
> >> x1 = 2.219 + (- 5.53545630592774e-06*(6.1 - 6.1253) - ...
> 9.41932297959802e-06)/0.0258357008522912
>
> x1 =
> 2.21864083513785
>
> Now, let me try it out in myfun.
>
> >> myfun([x1,6.1],eta,Dmin,Dmax)
> ans =
> 4.24534185583525e-10
>
> Indeed, this is pretty close to zero.
>
> Have I ever touched fsolve? NO. Look at your
> problem. Think about it. Don't just blindly assume
> the computer will be smart enough to solve your
> problem with no thought of yours invested. This
> only works in the movies or on some mindless TV
> show.
>
> John

Tags for this Thread

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us