Thread Subject: Function from 0 to 1

Hey guys,

Again with the same question on finding a function from 0 to 1. I have been trying to get this as closest to a linear function but being an s-curve. Always increasing with a value of 0 at 0 and a value of 1 at 3000. So far the best I got is,

l = @(x) gammainc(x/325,5);
ezplot(l,[0,3000])

However at 3000 it is not 1. If anyone can please help me find a function that is a little more linear than this one, and that is 1 at 3000.

Thanks all!

"Diego Zegarra" <diegozbb@gmail.com> wrote in message <gp906p$h4j$1@fred.mathworks.com>...
> Hey guys,
>
> Again with the same question on finding a function from 0 to 1. I have been trying to get this as closest to a linear function but being an s-curve. Always increasing with a value of 0 at 0 and a value of 1 at 3000. So far the best I got is,
>
> l = @(x) gammainc(x/325,5);
> ezplot(l,[0,3000])
>
> However at 3000 it is not 1. If anyone can please help me find a function that is a little more linear than this one, and that is 1 at 3000.
>
> Thanks all!

Yes, but you wanted an asymptote. An asymptotic
function never reaches the asymptote.

How do you define a "little more linear"? Is that
like a little more pregnant? Please state your goal
clearly. For example, a simple set of goals might
be this:

1. f(0) = 0
2. f'(0) = 0
3. f(3000) = 1
4. f'(3000) = 0

Even that is insufficient to define the shape of
an interpolant fully, as you should appreciate
from the last time you asked this question. It is
however, sufficient to define a cubic polynomial.
The inflection point will occur midway in the
interval of interest, so at x = 1500.

Perhaps it is instructive to tech you how one
might derive a polynomial with the desired
properties. Best is to solve the dimensionless
but related problem:

1. g(0) = 0
2. g'(0) = 0
3. g(1) = 1
4. g'(1) = 0

where

g(u) = a*u^3 + b*u^2 + c*u + d

and

f(x) = g(x/3000)

See that when x = 0, x/3000 = 0, but when
x = 3000, x/3000 = 1.

If f(0) = f'(0) = 0, then we can quickly see
that c = d = 0. I'll leave this part to you to
prove. (Do so.) All that is required is to evaluate
the polynomial, g(u), at u = 0. Do you see that
this gives d = 0? Next, differentiate the poly,
and substitute 0 into that. What is c?

So we know that

g(u) = a*u^3 + b*u^2

The derivative is

g'(u) = 3*a*u^2 + 2*b*u

Substitute u = 1 into both of them, and solve
for a and b.

   1 = a*1 + b*1 = a + b

   0 = 3*a*1 + 2*b*1 = 3*a + 2*b

Solving, we get

  a = -2

and

  b = 3

Now, put it all together.

g(u) = -2*u^3 + 3*u^2
f(x) = g(x/3000)

Write it in matlab as a function.

fun = @(x) polyval([-2 3 0 0],x/3000);
ezplot(fun,[0 3000])

See that the point of inflection occurs at x = 1500
as I predicted. Unfortunately, the location of that
inflection point is not controllable, nor is the shape.

A question for the interested student is, can we do
the same thing with a 4th order polynomial? See
that this would leave us with one free parameter.
Could we control the location of the point of
inflection with that free parameter? For example,
define the location of the inflection point as p,
where p is known to lie in the interval [0,1].

Would the polynomial g(u) still have good properties
over the interval [0,1]? What properties would be
appropriate? What constraints must apply to the
parameter (p) for this to happen?

John

"John D'Errico" <woodchips@rochester.rr.com> wrote in message <gp99jb$12l$1@fred.mathworks.com>...

> A question for the interested student is, can we do
> the same thing with a 4th order polynomial? See
> that this would leave us with one free parameter.
> Could we control the location of the point of
> inflection with that free parameter? For example,
> define the location of the inflection point as p,
> where p is known to lie in the interval [0,1].
>
> Would the polynomial g(u) still have good properties
> over the interval [0,1]? What properties would be
> appropriate? What constraints must apply to the
> parameter (p) for this to happen?

As a spoiler, if you want to solve the above
problem, don't read what lies below.

John



m lies in the interval [1/3 , 2/3] for G to be a
monotonic function on the interval [0,1].

m = 1/3;
a = 3*(2*m - 1)/(6*m^2 - 6*m + 1)
b = -2 - 2*a
c = 1 - a - b

G = [a b c 0 0]
subplot(1,2,1)
ezplot(@(x) polyval(G,x/3000),[0 3000])

John,

I would really appreciate if you could give me a reference where I can get this stuff. I need to do some curve fitting and do not want to be asking all the time if I can do it myself. Your explanation for all this was great by the way, very detailed and easy to follow. I am trying to see how you came with those numbers in the 4th level polynomial though.

FYI. I am trying to fit curves with same characteristics but having the inflection point before, also having it after and also being not so linear, a little more 'pregnant'. If you could give me a hint where I can learn this I would really appreciate it.

Thanks!

And my question on the 4th order polynomial is where is m coming from, that the part I dont get.

Thanks

"Diego Zegarra" <diegozbb@gmail.com> wrote in message <grlbba$3a2$1@fred.mathworks.com>...
> John,
>
> I would really appreciate if you could give me a reference where I can get this stuff. I need to do some curve fitting and do not want to be asking all the time if I can do it myself. Your explanation for all this was great by the way, very detailed and easy to follow. I am trying to see how you came with those numbers in the 4th level polynomial though.
>
> FYI. I am trying to fit curves with same characteristics but having the inflection point before, also having it after and also being not so linear, a little more 'pregnant'. If you could give me a hint where I can learn this I would really appreciate it.
>
> Thanks!

Rats. So now I need to remember what I did
then. 8-) Luckily it is just some basic calculus
as I recall.

Ok, so your question is where did m come from,
and I assume, how did I know that m has limits
[1/3, 2/3] for monotonic G(x) on the interval
[0,1]?

Start by applying the simple constraints

G(0) = 0
G(1) = 1
G'(0) = G'(1) = 0

If we assume initially that G(x) is a general 4th
order polynomial,

G(x) = a*x^4 + b*x^3 + c*x^2 + d*x + e

Then we know immediately that d = e = 0.
Why? What is the Taylor series for G, expanded
around x = 0? Remember, we know that G(0)=0.
That tells us that e = 0. Likewise, G'(0) = 0
implies that d = 0 in the above polynomial.

See what happens when we substitute 1 into
G now.

  G(1) = a + b + c = 1

This gives us one linear constraint on the
parameters. Look at the derivative.

  G'(x) = 4*a*x^3 + 3*b*x^2 + 2*c*x

so

  G'(1) = 4*a + 3*b + 2*c = 0

Finally, let us choose the location of the
inflection point. I called the point of
inflection m. What is an inflection point?
It is a point where the second derivative
changes sign.

  G''(x) = 12*a*x^2 + 6*b*x + 2*c

Thus we have

  6*a*m^2 + 3*b*m + c = 0

For this to have a solution, what might we
know about the coefficients? Look at the
discriminant of the quadratic. For a real
root to exist for the quadratic, what can we
say?

  9*b^2 - 24*a*c >= 0

Put it all together now.

(1) a + b + c = 1
(2) 4*a + 3*b + 2*c = 0
(3) 6*a*m^2 + 3*b*m + c = 0

Solve for a, b, c, in terms of m, I get these
relations:

  c = 1 - a - b
  b = -2 - 2*a

and, finally,

  3 - 6*m + a - 6*a*m + 6*a*m^2 = 0

This gives us a, as a function of m.

  a = 3*(2*m-1)/(6*m^2 - 6*m + 1)

Given a, we can now back out the values for b
and then c.

Look back at the discriminant expression. I'll use
my sympoly toolbox.

sympoly a b c
subs(subs(9*b^2 - 24*a*c,'c',1-a-b),'b',-2-2*a)
ans =
    36 + 12*a^2

Obviously, this is positive for all real a. So as it
turns out, as long as we can find a valid solution
for a above in terms of m, it looks like things
should work. Recall the relation for a.

  a = 3*(2*m-1)/(6*m^2 - 6*m + 1)

This has singularities at two points. So m cannot
go outside of that interval.

roots([6 -6 1])
ans =
      0.78868
      0.21132

To get a feeling for what happens, try this:

m = 1/3;
a = 3*(2*m - 1)/(6*m^2 - 6*m + 1);
b = -2 - 2*a;
c = 1 - a - b;
G = [a b c 0 0];
ezplot(@(x) polyval(G,x),[0 1])

Look at the resulting shape of the curve. Try
it again for a different value of m, perhaps 2/3.
Now, try it one more time, at m = 0.3. This is
just a wee bit lower than 1/3. What happened?

In fact, there was a SECOND inflection point in
the curve. That second inflection point in the
interval forces the curve to be non-monotonic.

  G''(x) = 12*a*x^2 + 6*b*x + 2*c

At what point will that second inflection point
drop down below 1?

gpp = subs(subs(12*a*x^2 + 6*b*x + 2*c,'c',1-a-b),'b',-2-2*a)
gpp =
    6 - 12*x + 2*a - 12*a*x + 12*a*x^2

subs(gpp,'x',1)
ans =
    -6 + 2*a

So, G''(1) == 0 when 2*a - 6 == 0, or when
a = 1/3. For any smaller values of a, the second
inflection point on this curve creeps inside the
interval [0,1].

We can do the same analysis when m is near 2/3,
as the second inflection point will move above 0.

It is all just calculus and a little algebra. If you are
unsure what is happening, make some plots, think
about what you see. To be honest, there are surely
other ways to have derived these relations, and I
have no certainty that I used the same logic some
time ago when I first did the analysis. Luckily I
came to the same conclusion. 8-)

John

"John D'Errico" <woodchips@rochester.rr.com> wrote in message <grlqpp$csu$1@fred.mathworks.com>...

> subs(gpp,'x',1)
> ans =
> -6 + 2*a
>
> So, G''(1) == 0 when 2*a - 6 == 0, or when
> a = 1/3. For any smaller values of a, the second
> inflection point on this curve creeps inside the
> interval [0,1].

Oh. Sorry. I had to run out before.

a = 1/3 will imply a value for m, such that the
curve has two inflection points in the interval.
It must be at m = 1/3.

John

John,

Thank you so much for ALL your invaluable help!! It is all coming back now, all that calculus, concavity tests, critical and inflection points and so much more. I have gone over all this and analyzed it using a calc book and have understand all the steps you did. I thought I was seriously never going to use that calc again, I guess you proved me wrong.

Thanks again!