Find a function with certain requirements

Question

I'm trying to find a function $y=f(x)$ that can be described as follows:

$f(x) = g(x) + c/(x-x_a)$.

With $f(x)$ I want to design a function with the following properties:

• $f(0) = 0$;
• $f(x)$ has a maximum at $x_1, \quad x_1>0$;
• the maximum of $f(x)$ is $\max(f(x))=f(x_1)=y_\max$
• $f(x_2) = 0, \quad x_2>x_1$;
• $f(x)$ has a vertical asymptote at $x=x_a, \quad x_a<0$;
• $f(x)$ has a second asymptote which is a function $g(x)$;

Thus: I want to set $x_1, x_2$ and $y_\max$ and provide a function $g(x)$ to create a function with the listed requirements.

I have tried to implement this by solving a system of non-linear equations, but I could not achieve a function with all the listed conditions being met. I used MATLAB function fsolve to solve the unknown parameters of $g(x), c$ and $x_a$, by solving $f(x)$ for the unknown parameters manually (this could be done using the symbolic math toolbox as well).

An example, to give you a better feeling for what I'm looking for, is shown in the figure below. In the example, the asymptote function I used was $g(x) = ax+b$, so $a$ and $b$ are additional parameters to solve.

From the figure you can see that conditions $f(0)=0$ and $f(x_2)=0, x_2=7$ have been met, but in this case I didn't manage to set the desired value of the maximum ($y_\max$) at the desired $x$-value ($x_1$). How do I do this (by hand and/or using MATLAB)? Do I need to solve $\mathrm{d}f(x)/\mathrm{d}x=0$ for the parameters involving the peak height and $x$-location, or can this be done in a different or more efficient way?

Answer 1

OK, so my initial answer totally missed the point, but after having considered this a bit more, I came up with the following: $$ \begin{align} f(x)&=-\frac23x+\frac{50}9+\frac{-200/27}{x+4/3} \end{align} $$ which appears to solve the exact problem you posed. Here is a Wolfram|Alpha table confirming that $f(0)=f(7)=0$ and $f(2)=2$ as desired.

---

Here is how I found it:

The blue curve is given by $(x+k)^2(y+1)=(k+2)^2$ passing through $(2,0)$ and having asymptotes $x=-k$ and $y=-1$. The purple curve is given by $$ h(x)=-\frac{(k+2)^2}{x+k}-x $$ and is a primitive function of the blue curve $y=\dfrac{(k+2)^2}{(x+k)^2}-1$. Then writing out $h(0)=h(7)$ and solving for $k$ yields $k=4/3$. With this we are able to define the red curve, the solution to your problem, namely by scaling and transposing $h$ as follows: $$ f(x)=\frac{2}{h(2)-h(0)}\cdot(h(x)-h(0)) $$ in which $h(2)-h(0)=3$ and $h(0)=-25/3$, and the constants I wrote in the beginning of my answer follows as an immediate consequence of this after a couple of simple calculations.

Answer 2

In this case, you just multiply the answer you have by a constant. If we call your current answer $f$, then you have $f(x_1) = c$ for some value $c$. Yo'd like it to be, say, $e$. So define $$ F(x) = \frac{e}{c} f(x) $$ and you've got your function.

Answer 3

This is an interpolation problem with conditions on the first and second derivatives.

Conditions on $f$: $f(x_0) = 0$, $f(x_1) = y_{\mbox{max}}$, $f(x_2) = 0$, for $x_0 = 0 < x_1 < x_2$

Conditions on $f'$: $f'(x_1) = 0$, $f'(x_a) = m_a$, $f'(x_b) = m_b$, for $x_a < x_0 < x_1 < x_2 < x_b$

Conditions on $f''$: $f''(x_1) < 0$

An easy way to find such an $f$ is assuming a polynomial solution with 7 unknown coefficients for the 7 conditions: $$ p(x) = \sum_{k=0}^6 a_k x^k \\ p'(x) = \sum_{k=0}^5 (k+1) a_{k+1} x^k \\ p''(x) = \sum_{k=0}^4 (k+1)(k+2) a_{k+2} x^k $$

This gives a system $$ \left( \begin{matrix} p(0) \\ p(x_1) \\ p(x_2) \\ p'(x_1) \\ p'(x_a) \\ p'(x_b) \\ p''(x_1) \end{matrix} \right) = \left( \begin{matrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & x_1 & x_1^2 & x_1^3 & x_1^4 & x_1^5 & x_1^6 \\ 1 & x_2 & x_2^2 & x_2^3 & x_2^4 & x_2^5 & x_2^6 \\ 0 & 1 & 2 x_1 & 3 x_1^2 & 4 x_1^3 & 5 x_1^4 & 6 x_1^5 \\ 0 & 1 & 2 x_a & 3 x_a^2 & 4 x_a^3 & 5 x_a^4 & 6 x_a^5 \\ 0 & 1 & 2 x_b & 3 x_b^2 & 4 x_b^3 & 5 x_b^4 & 6 x_b^5 \\ 0 & 0 & 2 & 6 x_1 & 12 x_1^2 & 20 x_1^3 & 30 x_1^4 \end{matrix} \right) \left( \begin{matrix} a_0 \\ a_1 \\ a_2 \\ a_3 \\ a_4 \\ a_5 \\ a_6 \end{matrix} \right) = \left( \begin{matrix} 0 \\ y_{\mbox{max}} \\ 0 \\ 0 \\ m_a \\ m_b \\ n \end{matrix} \right) \iff \\ M a = b $$ It depends on the values of $x_a,x_1, x_2,x_b$ if $M$ is invertible or not.

Example:

Note: This example session uses the open source Maxima computer algebra system.

The matrix $M$:

(%i) M : matrix([1,0,0,0,0,0,0],
[1,x1,x1^2,x1^3,x1^4,x1^5,x1^6],
[1,x2,x2^2,x2^3,x2^4,x2^5,x2^6],
[0,1,2*x1,3*x1^2,4*x1^3,5*x1^4,6*x1^5],
[0,1,2*xa,3*xa^2,4*xa^3,5*xa^4,6*xa^5],
[0,1,2*xb,3*xb^2,4*xb^3,5*xb^4,6*xb^5],
[0,0,2,6*x1,12*x1^2,20*x1^3,30*x1^4]);

                [ 1  0    0      0      0       0       0    ]
                [                                            ]
                [          2      3      4       5       6   ]
                [ 1  x1  x1     x1     x1      x1      x1    ]
                [                                            ]
                [          2      3      4       5       6   ]
                [ 1  x2  x2     x2     x2      x2      x2    ]
                [                                            ]
                [                  2      3       4       5  ]
(%o)            [ 0  1   2 x1  3 x1   4 x1    5 x1    6 x1   ]
                [                                            ]
                [                  2      3       4       5  ]
                [ 0  1   2 xa  3 xa   4 xa    5 xa    6 xa   ]
                [                                            ]
                [                  2      3       4       5  ]
                [ 0  1   2 xb  3 xb   4 xb    5 xb    6 xb   ]
                [                                            ]
                [                          2       3       4 ]
                [ 0  0    2    6 x1   12 x1   20 x1   30 x1  ]

Example settings:

(%i) [xa, x1, x2, xb] : [-1,2,7,8];
(%o)                           [- 1, 2, 7, 8]

which lead to

(%i) M, numer;
                    [ 1  0   0    0    0      0      0    ]
                    [                                     ]
                    [ 1  2   4    8    16    32      64   ]
                    [                                     ]
                    [ 1  7  49   343  2401  16807  117649 ]
                    [                                     ]
(%o)                [ 0  1   4   12    32    80     192   ]
                    [                                     ]
                    [ 0  1  - 2   3   - 4     5     - 6   ]
                    [                                     ]
                    [ 0  1  16   192  2048  20480  196608 ]
                    [                                     ]
                    [ 0  0   2   12    48    160    480   ]

Determinant:

(%i) determinant(M),numer;
(%o)                            181993392000

which is non-zero so this $M$ is invertible.

Target vector $b$:

(%i) b : matrix([0],[ymax],[0],[0],[ma],[mb],[n]);
                                   [  0   ]
                                   [      ]
                                   [ ymax ]
                                   [      ]
                                   [  0   ]
                                   [      ]
(%o)                               [  0   ]
                                   [      ]
                                   [  ma  ]
                                   [      ]
                                   [  mb  ]
                                   [      ]
                                   [  n   ]

Example values:

(%i) [ymax,ma,mb,n] : [2,1,-1,-1];
(%o)                         [2, 1, - 1, - 1]
(%i) b, numer;
                                    [  0  ]
                                    [     ]
                                    [  2  ]
                                    [     ]
                                    [  0  ]
                                    [     ]
(%o)                                [  0  ]
                                    [     ]
                                    [  1  ]
                                    [     ]
                                    [ - 1 ]
                                    [     ]
                                    [ - 1 ]

Solution vector $x$:

(%i) a : invert(M) . b, numer;
                           [          0.0          ]
                           [                       ]
                           [   1.739954605121048   ]
                           [                       ]
                           [  - 0.085881793862054  ]
                           [                       ]
(%o)                       [  - 0.23119708223252   ]
                           [                       ]
                           [   0.050675057740558   ]
                           [                       ]
                           [ - 0.0030721632354652  ]
                           [                       ]
                           [ 1.0983168004253304E-5 ]

This is the polynomial:

(%i) p(x) := a[1] + a[2]*x + a[3]*x^2 + a[4]*x^3 + a[5]*x^4 + 
a[6]*x^5 + a[7]*x^6;
                                 2       3       4       5       6
(%o)     p(x) := a  + a  x + a  x  + a  x  + a  x  + a  x  + a  x
                  1    2      3       4       5       6       7

First derivative

(%i) dp(x) := a[2] + 2*a[3]*x + 3*a[4]*x^2 + 4*a[5]*x^3 + 
5*a[6]*x^4 + 6*a[7]*x^5;
                                  2         3         4         5
(%o) dp(x) := a  + 2 a  x + 3 a  x  + 4 a  x  + 5 a  x  + 6 a  x
               2      3        4         5         6         7

Second derivative

(%i) d2p(x) := 2*a[3] + 6*a[4]*x + 12*a[5]*x^2 + 20*a[6]*x^3 +
30*a[7]*x^4;
                                      2          3          4
(%o) d2p(x) := 2 a  + 6 a  x + 12 a  x  + 20 a  x  + 30 a  x
                  3      4         5          6          7

With this we can check the solution visually:

(%i) plot2d(p(x),[x,-1,8]);
(%o) 

(%i) plot2d(dp(x),[x,-1,8]);
(%o)

(%i) plot2d(d2p(x),[x,-1,8]);
(%o)

The plots show that the conditions are satisfied by $p$.

Another nice method is interpolating piecewise with splines, see Spline interpolation.

Answer 4

Without loss of generality (by rescaling $x$ and $y$), one may assume $x_1=t$ (with $0<t<1$), $x_2=1$, and $y_{max}=1$.

Your parametric model is (setting $\frac{x_a}{x_2}=d$): $$f(x)=ax+b+\frac c{x-d},$$ which can be rewritten as $$ax^2+b'x+c'+df(x)=xf(x),$$ where $b'=b-ad,c'=c-bd$.

The derivative is such that

$$2ax+b'+df'(x)=f(x)+xf'(x).$$

From the first condition $f(0)=0$, you get $c'=0$. Plugging the remaining conditions, you get an easy linear system of $3$ equations in $3$ unknowns:

$$\begin{align}f(t)=1&\implies at^2+b't+d=t,\\ f'(t)=0&\implies 2at+b'=1,\\ f(1)=0&\implies a+b'=0.\end{align}$$

Then $$a=\frac1{2t-1},\\ b'=-\frac1{2t-1},\\ d=\frac{t^2}{2t-1}. $$ From there, $$b=b'+ad=\frac{(t-1)^2}{(2t-1)^2},\\ c=bd=\frac{t^2(t-1)^2}{(2t-1)^3}.$$

For $t<\frac12$ (i.e. $x_1<\frac{x_2}2$ like in your example), the value of $d$ (and $x_a$) is negative and this model works.