MATLAB Central Newsreader - solution of a transcendental equation

solution of a transcendental equation

Fulya — Mon, 14 Jan 2008 20:29:14 -0500

Hi all,

I have a transcendental equation:

y=a.*((b+c.*y.*x)./(b+c.*x))^(1-1./c)

I have to take y to the left side of equation and obtain an expression y=f(x). x is an array of 50 elements and a,b and c are constants.
Do you have any idea? What can I do?
b and c parameters are constantly changing.
I tried 'solve' command but it didn't do like I want.
I know that c>1.
Thanks a lot.

Fulya

Re: solution of a transcendental equation

Roger Stafford — Tue, 15 Jan 2008 00:25:18 -0500

Fulya <fulya_bagci@yahoo.com> wrote in message
<9544552.1200342600658.JavaMail.jakarta@nitrogen.mathforum.org>...
> Hi all,
>
> I have a transcendental equation:
>
> y=a.*((b+c.*y.*x)./(b+c.*x))^(1-1./c)
>
> I have to take y to the left side of equation and obtain an expression y=f(x).
x is an array of 50 elements and a,b and c are constants.
> Do you have any idea? What can I do?
> b and c parameters are constantly changing.
> I tried 'solve' command but it didn't do like I want.
> I know that c>1.
> Thanks a lot.
>
> Fulya
---------
  Possibly you can use matlab's 'fzero' function if you can come up with a
scheme for making the x0 estimates that function requires. If your b and c
values change by small amounts, you can probably use the solutions to the
previous b and c values for estimating the next b and c solutions.

  I don't have 'fzero' on my ancient system but I think it only works with
scalar-valued functions, so you may have to resort to a 50-step for-loop to
treat each value in the x vector separately in determining the corresponding y
vector, for given a, b, and c values.

  I notice that you can easily solve for x in the equation as a function of y,
which means there is only one value of x for each y-value (for a given a, b, c
set.) However, the reverse does not appear to be true. Given values of x may
correspond to more than one y value. In this case the x0 estimate for 'fzero'
will have to be carefully chosen so as to lead to the correct y. Either that or
you need to invoke some constraints on the y-value range (using two-
element x0's.)

  The fact that x can be expressed as a function of y makes it possible for you
to do some initial plotting for various values of your parameters a, b, and c,
so as to help with the proper choice of the estimator x0 values.

Roger Stafford

Re: solution of a transcendental equation

John D'Errico — Tue, 15 Jan 2008 01:26:02 -0500

Fulya <fulya_bagci@yahoo.com> wrote in message
<9544552.1200342600658.JavaMail.jakarta@nitrogen.mathforum.org>...
> Hi all,
>
> I have a transcendental equation:
>
> y=a.*((b+c.*y.*x)./(b+c.*x))^(1-1./c)
>
> I have to take y to the left side of equation and obtain an expression y=f(x).
x is an array of 50 elements and a,b and c are constants.
> Do you have any idea? What can I do?
> b and c parameters are constantly changing.
> I tried 'solve' command but it didn't do like I want.
> I know that c>1.
> Thanks a lot.

% Pick some arbitrary numbers for a, b, c.
a = 2;
b = 3;
c = pi;
% and a few random numbers for x.
n = 10;
x = rand(n,1);

% first, do it in a loop
yi = nan(size(x));
for i = 1:n
  funi = @(y) y-a.*abs((b+c.*y.*x(i))./(b+c.*x(i))).^(1-1./c);
  yi(i) = fzero(funi,1);
end

% verify the solutions, using a vectorized form
fun = @(y) y-a.*abs((b+c.*y.*x)./(b+c.*x)).^(1-1./c);

% solve all at once using fsolve
opts = optimset('fsolve');
opts.JacobPattern = speye(n,n);
yvec = fsolve(fun,ones(n,1),opts);

[yi,yvec,yi-yvec,fun(yi)]
ans =
       3.4128 3.4128 -5.3715e-08 0
       3.1859 3.1859 -1.6924e-09 0
       2.7411 2.7411 -3.9524e-14 0
       2.6615 2.6615 8.8818e-16 0
       2.6786 2.6786 -2.6645e-15 4.4409e-16
       2.4596 2.4596 8.8818e-16 4.4409e-16
        2.775 2.775 -1.2124e-13 0
       2.7785 2.7785 -1.3634e-13 0
       3.2486 3.2486 -4.7774e-09 4.4409e-16
       3.2146 3.2146 -2.7467e-09 4.4409e-16

HTH,
John

Re: solution of a transcendental equation

Fulya — Wed, 16 Jan 2008 17:28:57 -0500

OK. I used 'fsolve'.

I do least squares fitting of this function.It gives such an error and although the fit is good when I enter good initial parameters, it is not smooth and it gives what I write as the initial parameters.Error:

Conditioning of Gradient Poor - Switching To LM method
Optimization terminated: search direction less than TolX.

my code is
.
funi = @(y) mean(abs(y-a.*abs(Fsat+K.*y.*x)./ (Fsat+K.*x)).^(1-1./K));
%yd is experimental y, yi is the theoretical y
yi = fsolve(funi,yd);
.
val = yd - yi

second window:
read experimental x and yd data
[solution]=lsqnonlin(@OAOLmyLu,[8.0,3.0])
plot(x,yd)
plot(x,yi)

I thank individually to whow writes to me.

Re: solution of a transcendental equation

John D'Errico — Wed, 16 Jan 2008 18:59:01 -0500

Fulya <fulya_bagci@yahoo.com> wrote in message
<22495104.1200504568448.JavaMail.jakarta@nitrogen.mathforum.org>...
> OK. I used 'fsolve'.

Fsolve is inappropriate to solve a least squares
problem as you are doing.

> I do least squares fitting of this function.It gives such an error and although
the fit is good when I enter good initial parameters, it is not smooth and it
gives what I write as the initial parameters.Error:
>
> Conditioning of Gradient Poor - Switching To LM method

This is not an error message. It is a statement
about the algorithm's progress.

> Optimization terminated: search direction less than TolX.

This is not an error message. It is a statement
that the algorithm has finished its work.

> my code is
> .
> funi = @(y) mean(abs(y-a.*abs(Fsat+K.*y.*x)./ (Fsat+K.*x)).^(1-1./K));
> %yd is experimental y, yi is the theoretical y
> yi = fsolve(funi,yd);
> .
> val = yd - yi

What are you solving for? y? Or are you
trying to solve for the parameters Fsat, K, a?

John

Re: solution of a transcendental equation

Roger Stafford — Wed, 16 Jan 2008 20:58:02 -0500

Fulya <fulya_bagci@yahoo.com> wrote in message
<22495104.1200504568448.JavaMail.jakarta@nitrogen.mathforum.org>...
> OK. I used 'fsolve'.
>
> I do least squares fitting of this function.It gives such an error and although
the fit is good when I enter good initial parameters, it is not smooth and it
gives what I write as the initial parameters.Error:
>
> Conditioning of Gradient Poor - Switching To LM method
> Optimization terminated: search direction less than TolX.
>
> my code is
> .
> funi = @(y) mean(abs(y-a.*abs(Fsat+K.*y.*x)./ (Fsat+K.*x)).^(1-1./K));
> %yd is experimental y, yi is the theoretical y
> yi = fsolve(funi,yd);
> .
> val = yd - yi
>
> second window:
> read experimental x and yd data
> [solution]=lsqnonlin(@OAOLmyLu,[8.0,3.0])
> plot(x,yd)
> plot(x,yi)
>
> I thank individually to whow writes to me.
----------
  Fulya, I don't understand how your function 'funi' in this second article
relates to the original implicit function you defined in your first article:

y=a*((b+c*y*x)/(b+c*x))^(1-1/c)

Presumably 'Fsat' refers to 'b' and 'K' to 'c'. However, in this second article the
exponent, (1-1/c), operates on the difference between the two sides of the
equation, rather than simply to the fraction (b+c*y*x)/(b+c*x). Your 'funi' is
also not the same as John's 'funi' or his 'fun'. How do you account for the
differences here?

  Also you don't state anywhere whether quantity 'a' is less than or greater
than 1. This makes a very important difference to the appearance of the x-
versus-y curve. Also you don't state whether you are interested in the
portion where y > a or that where y < a. These are two separate pieces of the
curve. What kinds of values are you using for a, b, and c?

  As I mentioned earlier, x can be expressed as a function of y:

x = b/c*(t-1)/(y-t)

where

t = (y/a)^(c/(c-1)).

(Very likely you intend for y and a to both be positive quantities so as to get a
real value for t here.) You should try plotting this function from y ranging
from 0 up to something well above y = a to get a feeling for what you are
dealing with. You will be able to see the difference between the cases where
a > 1 and a < 1, as well as noting the asymptotes at y = 0 and y = a.

  If I were doing this problem (after doing some plotting) I would first use
'fzero' to solve it, before attempting to use 'fsolve'. Though 'fzero' can only
deal with scalar solutions one at a time, it is much simpler in concept than
'fsolve' and would allow you to better understand the problem you are dealing
with. After that you could try for better run time efficiency with 'fsolve'. You
should note that the results John got with 'fzero' gave solutions that appear
to be as accurate as matlab's double is capable of - that is, errors only out at
the least bit.

Roger Stafford

Re: solution of a transcendental equation

Fulya — Thu, 17 Jan 2008 08:07:21 -0500

I want to solve Fsat and K (refering to b and c respectively).a is a constant and 1.289 for this fit.The fitting equation is a trascendental equation so I must solve it first by 'fzero' or 'fsolve' before fitting with 'lsqnonlin'. (1-1/c) is the exponent of (b+c*y*x)/(b+c*x).y decreases from 1 to 0.4 . There are 180 experimental y values.I do least squares fitting by obtaining theoretical y values by fsolve for experimental x values. Fsat (b) and K(c)can be 2 and 6 respectively (reasonable values)
For a=1.289, b=2 and c=6 the curve starts from 1, goes flat a little and bends down to y=0.4 exponentially( it is a semilogx graphic).
I used fsolve rather than fzero because
yi = fsolve(funi,yexperimental)
There are 180 experimental y's, if I use fzero I must enter one(or two) parameter as x0 (fzero(funi,x0))(it does not accept more) and since theoretical y's must be close to experimental y's I chose fsolve command. x (experimental) starts from 0.004 and ends with 1.49.
It gives a statement as I wrote before and the fit is not smooth and it fits to what I write as the initial parameters for Fsat and K.

function val = OAOLmyLu(s)
global x yexperimental ytheoretical val
a=1.289;
funi = @(y) mean(abs(y-a.*abs((Fsat+K.*y.*x)./(Fsat+K.*x)).^(1-1./K)));
ytheoretical = fsolve(funi,yexperimental)
val = ytheoretical-yexperimental;

**
close all, clc
global x yexperimental ytheoretical val
read experimental x and y from txt file.
[solution]=lsqnonlin(@OAOLmyLu,[8.0,3.0]);
K=solution(1)
Fsat=solution(2)
calculate sum(val.^2)
semilogx(x, yexperimental)
semilogx(x, ytheoretical)

Re: solution of a transcendental equation

Roger Stafford — Thu, 17 Jan 2008 21:20:19 -0500

Fulya <fulya_bagci@yahoo.com> wrote in message
<5287596.1200557272380.JavaMail.jakarta@nitrogen.mathforum.org>...
> I want to solve Fsat and K (refering to b and c respectively).a is a constant
and 1.289 for this fit.The fitting equation is a trascendental equation so I
must solve it first by 'fzero' or 'fsolve' before fitting with 'lsqnonlin'. (1-1/c) is
the exponent of (b+c*y*x)/(b+c*x).y decreases from 1 to 0.4 . There are 180
experimental y values.I do least squares fitting by obtaining theoretical y
values by fsolve for experimental x values. Fsat (b) and K(c)can be 2 and 6
respectively (reasonable values)
> For a=1.289, b=2 and c=6 the curve starts from 1, goes flat a little and
bends down to y=0.4 exponentially( it is a semilogx graphic).
> I used fsolve rather than fzero because
> yi = fsolve(funi,yexperimental)
> There are 180 experimental y's, if I use fzero I must enter one(or two)
parameter as x0 (fzero(funi,x0))(it does not accept more) and since
theoretical y's must be close to experimental y's I chose fsolve command. x
(experimental) starts from 0.004 and ends with 1.49.
> It gives a statement as I wrote before and the fit is not smooth and it fits to
what I write as the initial parameters for Fsat and K.
>
> function val = OAOLmyLu(s)
> global x yexperimental ytheoretical val
> a=1.289;
> funi = @(y) mean(abs(y-a.*abs((Fsat+K.*y.*x)./(Fsat+K.*x)).^(1-1./K)));
> ytheoretical = fsolve(funi,yexperimental)
> val = ytheoretical-yexperimental;
>
> **
> close all, clc
> global x yexperimental ytheoretical val
> read experimental x and y from txt file.
> [solution]=lsqnonlin(@OAOLmyLu,[8.0,3.0]);
> K=solution(1)
> Fsat=solution(2)
> calculate sum(val.^2)
> semilogx(x, yexperimental)
> semilogx(x, ytheoretical)
--------
  Based on the data you have given us, Fulya, it looks as though you are
facing an impossible task, fitting the theoretical to the experimental data.
You stated that "x (experimental) starts from 0.004 and ends with 1.49." You
also stated that a > 1, that presumably b was to be positive, and that c was to
be greater than 1. Finally you stated that, "y decreases from 1 to 0.4 ".

  The code below uses a = 1.289 and the "reasonable values" b = 2 and c = 6
you suggested, and plots x versus y over the range .01 <= y <= a^c-.1 . It
uses the formula I mentioned earlier for finding x as a function of y, namely, t
= (y/a)^(c/(c-1)), x = b/c*(t-1)/(y-t), and makes a plot of the curve in yellow.
(You can easily check that x and y satisfy your original equation.) It also
includes the asymptotic lines y = 0 and y = a^c in red. You will note that
throughout this y range, x increases as y increases and the two lines y = 0
and y = a^c are approached asymptotically as x approaches minus and plus
infinity, respectively. (This is because the denominator, y-t, becomes zero
when y = 0 and when y = a^c.) If you use smaller limits on y, y=linspace
(1.293,2.547,n), so as to restrict x approximately to the experimental range
in question, you can see that these y values are very far away from the above-
mentioned experimental y values which range from 1 down to .4 and they are
moving in the wrong direction. If we use the latter experimental y-range, we
get negative values for theoretical x, which is again far away from the
experimental values.

  Assuming only that a > 1, b > 0, and c > 1, this formula of x as a function
of y shows that on any theoretical curve, x will be zero when t = 1 and
therefore when y = a. For larger values of x the y values will be larger.
Hence, for any such a, b, and c, the theoretical curve will always remain far
removed from the stated experimental y values and will be sloping in the
opposite direction.

  If all this is correct, you would have no chance of fitting a theoretical curve
of the kind you have defined to the experimental data, regardless of what
method you use to solve the equation.

% Calculate the theoretical curve above y = 0 and below y = a^c.
clear
a =1.289; b=2; c=6;
n = 1000;
y = linspace(.01,a^c-.1,n).';
t = (y/a).^(c/(c-1));
x = b/c*(t-1)./(y-t);
plot([x(1),x(n)],[a^c,a^c],'r-',[x(1),x(n)],[0,0],'r-',x,y,'y.')

Roger Stafford