# How to solve a differential equation when the right-hand side of the equation is a non-parametric curve?

I just asked another question (see here: How to calculate the shape of a curve given y coordinates and slope?) and was advised by the user who answered my question to ask a new question. I would appreciate any help on this topic. As in the previous thread, I apologise in advance if I'm not using the proper terminology, but I'm learning.

Although my problem is stated in the linked thread, I'll re-state it here with the details that required me to open a new thread.

My Question:

I would like to figure out the shape of a curve given the information in the following graph:

On the y axis, I'm showing the slope of the curve whose shape I'm trying to find. On the x axis, I'm showing the y-coordinate of the curve corresponding to that slope. I am missing information about the x-coordinates of my curve. The different dots are different measurements I have made in an experiment.

First, I fit a non-parametric curve to my data (in this case, a loess regression curve):

This gives me a non-parametric description of the relationship between dy/dx of the curve whose shape I'm trying to find out, and its y coordinates.

In this case, I can intuitively understand that my mystery curve will have a sigmoidal shape because at low and high values of y, the slope is small, and at intermediate values of y, the slope is high.

I just learnt, in the previous question, that my problem involves solving a differential equation. However, there are times when I don't have an equation describing the relationship between dx/dy and y (as in the previous question that I asked), but instead I have a non-parametric curve like a loess (local regression) curve or a spline.

How could I solve this problem?

• Your curve looks close to $y'=y(20-y)/10$ which can be solved as Bernoulli equation using $y=1/u$ so that $-u'=y'/y^2=2u-1/10$ resulting in $u(x)=1/20+(u_0-1/20)e^{-2x}$ and $$y(x)=\frac{20y_0}{y_0+(20-y_0)e^{-2x}}.$$ Is that about what you expected to find? – Lutz Lehmann May 29 '19 at 10:02
• To be fair, this was just an example dataset. My y' = f(y) plot could take any shape and I was looking for a general method to solve this problem rather than an analytical solution, which is why I liked JJacquelin's second method, which I think is valid? – Ender May 29 '19 at 10:51
• You can use something like the second approach. Look up your numerics package for interpolation and do something like f = interpolatorfunction(y_array, dydx_array). Then call a ODE integrator with this function f and some sensible initial value for your overall situation to get a table or piecewise polynomial approximation for y(x). – Lutz Lehmann May 29 '19 at 11:10

UPDATED ANSWER. I didn't delete my first answer (even not corresponding to the question) because this first answer gives a useful example to understand the general approach to solve a differential equation when some functions are not explicit but given on the form of data.

FIRST MEHOD :

Fit a polynomial equation with the data. The curve looks like a parabola. $$y'=\frac{dy}{dx}\simeq ay^2+by+c \tag 1$$ Usual linear regression leads to $$\begin{cases} a\simeq -0.096 \\ b\simeq 1.91 \\ c\simeq 1.26 \end{cases}$$

Note that the data is not accurate because it comes from a numerical scan of the graph given by Ender instead of a numerical table.

Of course a better fit should be obtained with a polynomial equation of higher degree. But this is of no interest with such a data of bad accuracy.

Integrating Eq.$$(1)$$ leads to the approximate solution of the differential equation : $$x(y)=\int\frac{dy}{y^2+by+c}+\text{constant} \tag 2$$ No boundary condition is specified in the question. As a consequence the integration constant cannot be determined.

However we need a fully determined function in order to proceed to further numerical calculus. Supposing that the condition is $$x(y_1)=0$$ the above solution would be : $$x(y)=\int_{y_1}^y\frac{d\zeta}{\zeta^2+b\zeta+c} \tag 3$$ Of course, if another boundary condition is specified, one have change it in the numerical calculus.

The result of the first method is represented in the next figure :

There is no need for the analytic expression of the integral $$(3)$$. Numerical integration is sufficient with the usual incremental method : $$x(y+\delta y)=x(y)+\frac{1}{ay^2+by+c}\delta y$$ The computation and drawing were carried out with $$\delta y=0.01$$ for example.

Nevertheless one can use the analytical solution : $$x(y)= \frac{2}{\sqrt{4ac-b^2}}\left(\tan^{-1}\left(\frac{2ay+b}{\sqrt{4ac-b^2}} \right) -\tan^{-1}\left(\frac{2ay_1+b}{\sqrt{4ac-b^2}} \right)\right)$$ Or the inverse function : $$y(x)=-\frac{b}{2a}\pm\frac{\sqrt{4ac-b^2}}{2a}\tan\left(\frac12\sqrt{4ac-b^2}\: x+\tan^{-1}\left( \frac{2ay_1+b}{\sqrt{4ac-b^2}}\right) \right)$$ This is a complicated formula. That is why, in practice, it is easier tu proceed with numerical integration as pointed out above. 

SECOND METHOD :

Numerical integration of the data $$(y_1,y'_1),\:...\:, (y_k,y'_k),\:...\:, (y_{20},y'_{20})$$

With the same assumption $$x(y_1)=0$$ than in the first method, again one would have to add a convenient constant if the boundary condition was different.

$$x_1=0\quad;\quad x_k=x_{k-1}+\frac12\left(\frac{1}{y'_{k-1}}+\frac{1}{y'_k}\right)(y_k-y_{k-1})\qquad [\:2\leq k\leq 20\:]$$

On the next figure the result of the second method (in blue) is compared to the result of the first method (in red).

The discrepancy is probably due to the approximative method of integration in the second method. They are not enough points which are not regularly distributed. Some of the gaps between two consecutive points are too large.

That is why the first method is recommended especially in case of small number of points.

Moreover the first method provides the result on a continuous form while the second method provides the result on discret form. Often it is advantageous for further calculus to have a continuous function instead of a data table.

• thanks for the update, JJacquelin, I'm setting this as the correct answer although I agree it's good to have the other one around for reference. – Ender May 30 '19 at 9:25
• Hi @Jjacquelin I was going through this again and I had a question - why are the calculations for the second method/trapezoidal approximation the same in both answers (just changing the y symbol for x in the answer below), but the numbers/answers are different? – Ender Jul 23 '19 at 15:38
• IN THE FIRST METHOD, a least mean squares regression is done first. Regression is not an exact process and some small deviations appear. Then the integration is not numeric but analytic. So there is no additional deviation due to a numerical integration which doesn't exist. THE SECOND METHOD is a numerical integration. The numerical integration isn't an exact process. The small deviations come from the numerical integration. COMPARING BOTH METHODS : Since the causes of the small deviations are not the same in each method one cannot expect that the numerical results be exactly identical. – JJacquelin Jul 23 '19 at 20:59
• Thanks. But I mean comparing the second method in this answer and the second method in your other answer on this page (see below), where you solved this problem for y' = f(x) instead of y' = f(y). The trapezoid approximation you used is exactly the same in both cases. It's just the symbol "x" or "y" that changes. However, your numerical answers to both are very different – Ender Jul 23 '19 at 21:10
• mathworld.wolfram.com/NumericalIntegration.html . In the second method the numerical integration is for the function $1/f(y)$ with respect to the variable $y$ $$x=\int \frac{1}{f(y)}dy$$ because $f(y)=\frac{dy}{dx}\quad\to\quad dx=\frac{1}{f(y)}dy$ . – JJacquelin Jul 24 '19 at 6:08

FIRST MEHOD :

Fit a polynomial equation with the data. The curve looks like a parabola. $$y'=\frac{dy}{dx}\simeq ax^2+bx+c$$ Usual linear regression leads to $$\begin{cases} a\simeq -0.096 \\ b\simeq 1.91 \\ c\simeq 1.26 \end{cases}$$

Note that the data is not accurate because it comes from a numerical scan of the graph given by Ender instead of a numerical table.

Of course a better fit should be obtained with a polynomial equation of higher degree. But this is of no interest with such a data of bad accuracy.

Integrating Eq.$$(1)$$ leads to the approximate solution of the differential equation : $$y=\frac{a}{3}x^3+\frac{b}{2}x^2+cx+C$$ No boundary condition is specified in the question. As a consequence the integration constant $$C$$ cannot be determined.

However we need a numerical value for $$C$$ in order to proceed to further numerical calculus. Supposing that the condition is $$y(x_1)=0$$ the constant would be $$C=-(\frac{a}{3}x_1^3+\frac{b}{2}x_1^2+cx_1)\simeq -0.258$$

Of course, if another boundary condition is specified, one have to compute the corresponding value of $$C$$ and change it in the next numerical calculus.

The result of the first method is represented in the next figure :

This result and the approximate equation $$y=\frac{a}{3}x^3+\frac{b}{2}x^2+cx+C$$ are valid on the range of $$x$$ determined by the original data, about $$0. Outside this range one cannot guaranty the correctness.

SECOND METHOD :

Numerical integration of the data $$(x_1,y'_1),\:...\:, (x_k,y'_k),\:...\:, (x_{20},y'_{20})$$

With the same assumption $$y(x_1)=0$$ than in the first method, again one would have to add a convenient constant if the boundary condition was different.

$$y_1=0\quad;\quad y_k=y_{k-1}+\frac12\left(y'_{k-1}+y'_k\right)(x_k-x_{k-1})\qquad [\:2\leq k\leq 20\:]$$

The results of both methods are very close. The curves are almost indistinguishable : drawn in red for the first method and blue for the second.

• At the point of your answer, was the graph labeled differently? As the horizontal axis is now labeled "$y$", the represented curve is $y'=f(y)$, so you would have to solve the (separable) Riccati equation $y'(x)=ay(x)^2+by(x)+c$. – Lutz Lehmann May 29 '19 at 9:53
• OK. My answer is for $y'=f(x)$, not for $y'=f(y)$ .I didn't noticed that on the Ender's graph. – JJacquelin May 29 '19 at 10:17
• Wouldn't the second method still be valid for y' = f(y)? – Ender May 29 '19 at 10:52
• No, not directly, as this second method is also for the case that the graph is of a function of $x$. But you could still use a piecewise linear interpolation of the data as an approximation of $f$ and then use any ODE solver, or chain the exponential functions that are solution of the linear pieces together to a continuous function. – Lutz Lehmann May 29 '19 at 11:00
• Yes, sure the method will work. But not as it is written in my answer. I will edit an updated answer as soon as I will have available time. – JJacquelin May 29 '19 at 11:00