"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')(1gammainc(xx(2)*Dmin,xx(1)+4, 'upper'))))^2./((gamma(xx(1)+5)*(gammainc(xx(2)*Dmax,xx(1)+5,'lower')(1gammainc(xx(2)*Dmin,xx(1)+5, 'upper'))))*(gamma(xx(1)+3)*(gammainc(xx(2)*Dmax,xx(1)+3,'lower')(1gammainc(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 mfile,
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.4193e06
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.53545630592774e06
>> myfun(xx0,eta,Dmin,Dmax)
ans =
9.41932297959802e06
So a VERY good LINEAR model for this function is
z = 9.41932297959802e06 + 0.0258357008522912*(x1  2.219) + 5.53545630592774e06*(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.53545630592774e06*(6.1  6.1253)  ...
9.41932297959802e06)/0.0258357008522912
x1 =
2.21864083513785
Now, let me try it out in myfun.
>> myfun([x1,6.1],eta,Dmin,Dmax)
ans =
4.24534185583525e10
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
