System of delay differential equations: using first interpolation as second initial condition

I am trying to solve numerically the following system of two coupled delay differential equations:

$ $ \dot x(t)=-\gamma x(t)-\frac{\gamma}{4}e^{i\omega_0\tau_1}y(t-\tau_1)\theta(t-\tau_1)+\frac{\gamma}{4}e^{i\omega_0\tau_2}y(t-\tau_2)\theta(t-\tau_2)+\frac{\gamma}{2}e^{i\omega_0\tau_3}x(t-\tau_3)\theta(t-\tau_3),$ $ $ $ \dot y(t)= -\frac{\gamma}{2}y(t)-\frac{\gamma}{4}e^{i\omega_0\tau_1}x(t-\tau_1)\theta(t-\tau_1)+\frac{\gamma}{4}e^{i\omega_0\tau_2}x(t-\tau_2)\theta(t-\tau_2).$ $ where $ \tau_1<\tau_2<\tau_3$ . The parameters $ \gamma, \omega_0$ are constants, and $ \theta(t)$ is the Heaviside step function. The history of the system is known for $ 0\leq t\leq\tau_1$ : $ $ x(t)=e^{-\gamma t}, y(t)=e^{-\gamma t/2}.$ $ Here what I tried:

I first solved the system for $ 0\leq t\leq\tau_2$ using the aforementioned initial history using NDSolve:

\[Gamma] = 1.0; \[Omega]0 = 2 Pi; \[Tau]1 = 1.0; \[Tau]2 = 2.0; \[Tau]3 = 3.0;  sol1 = NDSolve[{x'[   t] == - \[Gamma] x[t] - (\[Gamma]/4) E^(I \[Tau]1 \[Omega]0)     y[t - \[Tau]1],  y'[t] == - 0.5 \[Gamma] y[t] - (\[Gamma]/4) E^(    I \[Tau]1 \[Omega]0) x[t - \[Tau]1],  x[t /; t <= \[Tau]1] == (1.0/Sqrt[2.0]) Exp[-\[Gamma] t],  y[t /; t <= \[Tau]1] == (1.0/Sqrt[2.0]) Exp[-0.5 \[Gamma] t]}, {x,  y}, {t, 0, \[Tau]2}]; 

I get the following solution for $ |x(t)|^2$ and $ |y(t)|^2$ : abs[x]^2,abs[y]^2

The problem arises when I use this first interpolated solution as the initial history to solve for the next interval of time:

sol2 = NDSolve[{x'[   t] == - \[Gamma] x[t] - (\[Gamma]/4) E^(I \[Tau]1 \[Omega]0)     y[t - \[Tau]1] + (\[Gamma]/4) E^(I \[Tau]2 \[Omega]0)     y[t - \[Tau]2],  y'[t] == - 0.5 \[Gamma] y[t] - (\[Gamma]/4) E^(    I \[Tau]1 \[Omega]0) x[t - \[Tau]1] + (\[Gamma]/4) E^(    I \[Tau]2 \[Omega]0) x[t - \[Tau]2],  x[t /; t <= \[Tau]2] == Evaluate[x[t] /. sol1],  y[t /; t <= \[Tau]2] == Evaluate[y[t] /. sol1]}, {x, y}, {t,  0, \[Tau]3}];  

This time I get the following messages: enter image description here enter image description here

It seems that the second NDSolve (sol2) does not allow the interpolation of the first result as initial history. Any suggestion? Thank you in advance.

Interpolation of render with fixed physics update

In Fix Your Timestep, they briefly address the remainder time in relation to rendering. Saying that remainder can be used to generate an interpolated rendering.

We can use this remainder value to get a blending factor between the previous and current physics state simply by dividing by dt. This gives an alpha value in the range [0,1] which is used to perform a linear interpolation between the two physics states to get the current state to render.

http://web.archive.org/web/20190506030345/https://gafferongames.com/post/fix_your_timestep/

So let’s say your last two physics updates we’ll call A and B are from 00:10.033 and 00:10.066, and it is now 00:10.070 when we want to perform a render. We have a remainder of .004.

I take “interpolation between the two physics states” to mean we compare all the objects in update A and B and slide them back from B towards A by 88% (.033-.004)/.033). This would mean I’m actually rendering the physical state at 00:10.037, correct? So my physics updates A and B are really more like “previous” and “next” and my interpolation is between the previous and the next, correct?

Interpolation Inequality’s Proof

Let $ \Omega \subseteq R^{n}$ bounded domain. I need to prove that for $ u\in H^{2}(\Omega)\cap H^{1}_{0}(\Omega)$ : \begin{equation} \|\nabla u\|_{L^{2}(\Omega)}^{2}\leq \|u\|_{L^{2}(\Omega)}\|\Delta u\|_{L^{2}(\Omega)} \end{equation} I know that I should prove it for a $ u\in C^{2}_{0}(\Omega)$ and then use the Global Approximation Theorem with smooth functions to extend $ u\in C^{2}_{0}(\Omega)$ to $ u\in H^{2}(\Omega)\cap H^{1}_{0}(\Omega)$ . Also, I know the following expression could be useful: \begin{equation} \Delta(\frac{u^{2}}{2})=|\nabla u|^{2}+u\Delta u \end{equation}

Python multidimensional Lookup interpolation

To interpolate within a dataset, I want to create a multidimensional Lookup table in python. In this case, the tree input arrays x,y and z are given together with the output quantity a.

I was able to visualize the impact of x and y with a 2D heatmap:

xi = np.linspace(np.amin(x),np.amax(x),100) yi = np.linspace(np.amin(y),np.amax(y),100)  zi = griddata((x, y), a, (xi[None,:], yi[:,None]), method='linear')  CS = plt.contour(xi,yi,zi,15,linewidths=0.5,colors='k') CS = plt.contourf(xi,yi,zi,15,cmap=plt.cm.bwr_r) plt.colorbar() 

2D heatmap

However, my goal is not visualization. In the end I want to pass some x,y and z values to the Lookup and get an interpolated value of a.

I’ve allready found something related here, but the proposed function did not work. Is there any library I’ve missed so far or a suitable way to perform a multidimensional interpolation?

Lagrange interpolation formula


Give the formula of the $ 1$ st degree Lagrange polynomial $ L(x)$ interpolating a function $ f$ at the points $ 0$ and $ 1$ . Give the formula for the error $ L – f$ . Finally, show that

$ $ \sup_{x \in [0, 1]} |L(x) – f(x)| \leq \frac{1}{8} \sup_{[0, 1]} |f”|. $ $

To do this problem, I used the formula

$ $ L(x) = \sum_{i=0}^{n} l_{i}(x) y_{i},$ $

where $ $ l_{i}(x) = \prod_{0 \leq j \leq n \text{ and } j \neq i} \frac{x – x_{j}}{x_{i} – x_{j}}$ $

with $ x_{0} = 0, x_{1} = 1$ and $ y_{0} = f(x_{0}), y_{1} = f(x_{1})$ .

So first,

$ $ l_{0}(x) = \frac{x – x_{1}}{x_{0} – x_{1}} = \frac{x – 1}{-1} = 1 – x. $ $

Also,

$ $ \ell_{1}(x) = \frac{x – x_{0}}{x_{1} – x_{0}} = x .$ $

So I got $ L(x) = (1 – x)y_{0} + xy_{1}$ as my answer. Is this correct? I’m new to Lagrange interpolation and just wanna make I’m doing everything right.


Assuming this is right, I have the error formula provided in my book:

$ $ |L(x) – f(x)| \leq \frac{M_{n + 1}}{(n + 1)!} |\pi_{n + 1}(x)|,$ $

where $ M_{n + 1} = \max_{\xi \in [a, b]} |f^{(n + 1)}(\xi)|$ and $ \pi_{n + 1}(x) = \prod_{i=0}^{n} (x-x_{i})$ .

This looks really similar to what I need to do next, but I would like some guidance on how to do the next two parts.

Cubic spline interpolation in Python

I implemented the cubic spline interpolation explained in
https://en.wikipedia.org/wiki/Spline_interpolation as a Python class.

Following the best practices, I protected my code with tests (I preferred doctest for this example). In this way, we are good to refactor now.

Could you please help my to improve the quality of the code? You may just share your overall feedback or suggest some improvements.

Thanks,
Ivan

Code: CubicSplineStruct.py

""" Natural cubic spline interpolation Reference: https://en.wikipedia.org/wiki/Spline_interpolation To run doc tests: python -m doctest CubicSplineStruct.py  >>> cubicSplineStruct = CubicSplineStruct() >>> cubicSplineStruct.m_n = 4 >>> cubicSplineStruct.m_xvalues =  [0.0, 10./3., 20./3., 10.] >>> cubicSplineStruct.computeYtoKMatrix() >>> cubicSplineStruct.m_yvalues = [128., -64., 128., -64.] >>> cubicSplineStruct.computeKCoeffs() >>> print(cubicSplineStruct.interpolate(10./3.)) -64.0 >>> print(cubicSplineStruct.interpolate(5.)) 32.0 >>> print(cubicSplineStruct.interpolate(20./3.)) 128.0 """  import numpy as np from bisect import bisect_left  class CubicSplineStruct:     """     >>> cubicSplineStruct = CubicSplineStruct()     >>> cubicSplineStruct.m_n = 3     >>> cubicSplineStruct.m_xvalues =  [0., 10., 42.]     >>> cubicSplineStruct.computeYtoKMatrix()     >>> cubicSplineStruct.m_yvalues =   [100., 75., 95.]     >>> cubicSplineStruct.computeKCoeffs()     >>> y = cubicSplineStruct.interpolate(30.)     """     def __init__(self):         """         >>> cubicSplineStruct = CubicSplineStruct()         >>> print(cubicSplineStruct.m_n)         0         >>> print(cubicSplineStruct.m_xvalues)         []         >>> print(cubicSplineStruct.m_yvalues)         []         >>> print(cubicSplineStruct.m_kMatrix)         []         >>> print(cubicSplineStruct.m_yMatrix)         []         >>> print(cubicSplineStruct.m_ytoKMatrix)         []         >>> print(cubicSplineStruct.m_kCoeffs)         []         """         self.m_n = 0         self.m_xvalues = []         self.m_yvalues = []         self.m_kMatrix = np.matrix(np.zeros(shape=(0,0)))         self.m_yMatrix = np.matrix(np.zeros(shape=(0,0)))         self.m_ytoKMatrix = np.matrix(np.zeros(shape=(0,0)))         self.m_kCoeffs = []         pass      def pushFirstEquationToKMatrix(self, x0, x1):         """         >>> cubicSplineStruct = CubicSplineStruct()         >>> cubicSplineStruct.m_kMatrix = np.matrix(np.zeros(shape=(1,5)))         >>> cubicSplineStruct.pushFirstEquationToKMatrix(1.0, 1.5)         >>> print(cubicSplineStruct.m_kMatrix[0, 0]) # 2./(x1 - x0)         4.0         >>> print(cubicSplineStruct.m_kMatrix[0, 1]) # 1./(x1 - x0)         2.0         """         self.m_kMatrix[0, 0] = 2./(x1 - x0)         self.m_kMatrix[0, 1] = 1./(x1 - x0)      def pushLastEquationToKMatrix(self, xnm1, xn):         """         >>> cubicSplineStruct = CubicSplineStruct()         >>> cubicSplineStruct.m_kMatrix = np.matrix(np.zeros(shape=(5,5)))         >>> cubicSplineStruct.pushLastEquationToKMatrix(1.0, 1.5)         >>> print(cubicSplineStruct.m_kMatrix[-1, -1]) # 2./(xn - xnm1)         4.0         >>> print(cubicSplineStruct.m_kMatrix[-1, -2]) # 1./(xn - xnm1)         2.0         """         self.m_kMatrix[-1, -1] = 2./(xn - xnm1)         self.m_kMatrix[-1, -2] = 1./(xn - xnm1)      def pushMiddleEquationToKMatrix(self, i, xim1, xi, xip1):         """         >>> cubicSplineStruct = CubicSplineStruct()         >>> cubicSplineStruct.m_kMatrix = np.matrix(np.zeros(shape=(4,5)))         >>> cubicSplineStruct.pushMiddleEquationToKMatrix(3, 1.0, 1.5, 1.75)         >>> print(cubicSplineStruct.m_kMatrix[3, 2]) # 1./(xi - xim1)         2.0         >>> print(cubicSplineStruct.m_kMatrix[3, 3]) # 2./(xi - xim1) + 2./(xip1 - xi)         12.0         >>> print(cubicSplineStruct.m_kMatrix[3, 4]) # 1./(xip1 - xi)         4.0         """         self.m_kMatrix[i, i-1] = 1./(xi - xim1)         self.m_kMatrix[i, i] = 2./(xi - xim1) + 2./(xip1 - xi)         self.m_kMatrix[i, i + 1] = 1./(xip1 - xi)      def computeKMatrix(self):         """         >>> cubicSplineStruct = CubicSplineStruct()         >>> cubicSplineStruct.m_n = 4         >>> cubicSplineStruct.m_xvalues = [1.0, 1.5, 1.75, 2.25]         >>> cubicSplineStruct.computeKMatrix()         >>> print(cubicSplineStruct.m_kMatrix)         [[ 4.  2.  0.  0.]          [ 2. 12.  4.  0.]          [ 0.  4. 12.  2.]          [ 0.  0.  2.  4.]]         """         self.m_kMatrix = np.matrix(np.zeros(shape=(self.m_n, self.m_n)))         self.pushFirstEquationToKMatrix(self.m_xvalues[0], self.m_xvalues[1])         for i in range(1, self.m_n-1):             self.pushMiddleEquationToKMatrix(i, self.m_xvalues[i-1], self.m_xvalues[i], self.m_xvalues[i+1])         self.pushLastEquationToKMatrix(self.m_xvalues[-2], self.m_xvalues[-1])      def pushFirstEquationToYMatrix(self, x0, x1):         """         >>> cubicSplineStruct = CubicSplineStruct()         >>> cubicSplineStruct.m_yMatrix = np.matrix(np.zeros(shape=(1,5)))         >>> cubicSplineStruct.pushFirstEquationToYMatrix(1.0, 1.5)         >>> print(cubicSplineStruct.m_yMatrix[0, 0]) # -3./(x1 - x0)**2         -12.0         >>> print(cubicSplineStruct.m_yMatrix[0, 1]) # 3./(x1 - x0)**2         12.0         """         self.m_yMatrix[0, 0] = -3./(x1 - x0)**2         self.m_yMatrix[0, 1] = 3./(x1 - x0)**2      def pushLastEquationToYMatrix(self, xnm1, xn):         """         >>> cubicSplineStruct = CubicSplineStruct()         >>> cubicSplineStruct.m_yMatrix = np.matrix(np.zeros(shape=(5,5)))         >>> cubicSplineStruct.pushLastEquationToYMatrix(1.0, 1.5)         >>> print(cubicSplineStruct.m_yMatrix[-1, -1]) # 3./(xn - xnm1)**2         12.0         >>> print(cubicSplineStruct.m_yMatrix[-1, -2]) # -3./(xn - xnm1)**2         -12.0         """         self.m_yMatrix[-1, -1] = 3./(xn - xnm1)**2         self.m_yMatrix[-1, -2] = -3./(xn - xnm1)**2       def pushMiddleEquationToYMatrix(self, i, xim1, xi, xip1):         """         >>> cubicSplineStruct = CubicSplineStruct()         >>> cubicSplineStruct.m_yMatrix = np.matrix(np.zeros(shape=(4,5)))         >>> cubicSplineStruct.pushMiddleEquationToYMatrix(3, 1.0, 1.5, 1.75)         >>> print(cubicSplineStruct.m_yMatrix[3, 2]) # -3./(xi - xim1)**2         -12.0         >>> print(cubicSplineStruct.m_yMatrix[3, 3]) # 3./(xi - xim1)**2 - 3./(xip1 - xi)**2         -36.0         >>> print(cubicSplineStruct.m_yMatrix[3, 4]) # 3./(xip1 - xi)**2         48.0         """         self.m_yMatrix[i, i-1] = -3./(xi - xim1)**2         self.m_yMatrix[i, i] = 3./(xi - xim1)**2 - 3./(xip1 - xi)**2         self.m_yMatrix[i, i + 1] = 3./(xip1 - xi)**2      def computeYMatrix(self):         """         >>> cubicSplineStruct = CubicSplineStruct()         >>> cubicSplineStruct.m_n = 4         >>> cubicSplineStruct.m_xvalues = [1.0, 1.5, 1.75, 2.25]         >>> cubicSplineStruct.computeYMatrix()         >>> print(cubicSplineStruct.m_yMatrix)         [[-12.  12.   0.   0.]          [-12. -36.  48.   0.]          [  0. -48.  36.  12.]          [  0.   0. -12.  12.]]         """         self.m_yMatrix = np.matrix(np.zeros(shape=(self.m_n, self.m_n)))         self.pushFirstEquationToYMatrix(self.m_xvalues[0], self.m_xvalues[1])         for i in range(1, self.m_n-1):             self.pushMiddleEquationToYMatrix(i, self.m_xvalues[i-1], self.m_xvalues[i], self.m_xvalues[i+1])         self.pushLastEquationToYMatrix(self.m_xvalues[-2], self.m_xvalues[-1])      def computeYtoKMatrix(self):         """         Should be called when x knot values are updated          >>> cubicSplineStruct = CubicSplineStruct()         >>> cubicSplineStruct.m_n = 4         >>> cubicSplineStruct.m_xvalues = [0.0, 10./3., 20./3., 10.]         >>> cubicSplineStruct.computeYtoKMatrix()         >>> print(cubicSplineStruct.m_ytoKMatrix)         [[-0.38  0.48 -0.12  0.02]          [-0.14 -0.06  0.24 -0.04]          [ 0.04 -0.24  0.06  0.14]          [-0.02  0.12 -0.48  0.38]]         """         self.computeKMatrix()         self.computeYMatrix()         self.m_ytoKMatrix = self.m_kMatrix.I*self.m_yMatrix      def computeKCoeffs(self):         """         Should be called when y knot values are updated          >>> cubicSplineStruct = CubicSplineStruct()         >>> cubicSplineStruct.m_n = 2         >>> cubicSplineStruct.m_ytoKMatrix = np.mat('1 2; 4 5')         >>> cubicSplineStruct.m_yvalues =[3., 4.]         >>> cubicSplineStruct.computeKCoeffs()         >>> print(cubicSplineStruct.m_kCoeffs)         [11.0, 32.0]         """            kCoeffs = np.array(self.m_yvalues)*self.m_ytoKMatrix.T         self.m_kCoeffs = [kCoeffs[0, i] for i in range(self.m_n)]      def interpolateOnInterval(self, intervalIndex, x):         """         >>> cubicSplineStruct = CubicSplineStruct()         >>> cubicSplineStruct.m_xvalues = [None, 10., 42.]         >>> cubicSplineStruct.m_yvalues = [None, 128., -64.]         >>> cubicSplineStruct.m_kCoeffs = [None, 2., 6.]         >>> print(cubicSplineStruct.interpolateOnInterval(1, 10.))         128.0         >>> print(cubicSplineStruct.interpolateOnInterval(1, 18.))         98.0         >>> print(cubicSplineStruct.interpolateOnInterval(1, 34.))         -58.0         >>> print(cubicSplineStruct.interpolateOnInterval(1, 42.))         -64.0         """          x1 = self.m_xvalues[intervalIndex]         x2 = self.m_xvalues[intervalIndex+1]         y1 = self.m_yvalues[intervalIndex]         y2 = self.m_yvalues[intervalIndex+1]         t = (x - x1)/(x2 - x1)         a = computeACoeff(x1, x2, y1, y2, self.m_kCoeffs[intervalIndex])         b = computeBCoeff(x1, x2, y1, y2, self.m_kCoeffs[intervalIndex+1])         return (1-t)*y1 + t*y2 + t*(1-t)*(a*(1-t)+b*t)      def interpolate(self, x):         """         >>> cubicSplineStruct = CubicSplineStruct()         >>> cubicSplineStruct.m_xvalues = [None, 10., 42.]         >>> cubicSplineStruct.m_yvalues = [None, 128., -64.]         >>> cubicSplineStruct.m_kCoeffs = [None, 2., 6.]         >>> cubicSplineStruct.interpolate(18.)         98.0         """         if len(self.m_xvalues) == 0:             return 0.          intervalLowerBound = findLowerBound(self.m_xvalues, x)         return self.interpolateOnInterval(intervalLowerBound, x)  def computeACoeff(x1, x2, y1, y2, k):     """     >>> print(computeACoeff( 10., 42., 128., -64., 2.))     256.0     """     return k*(x2 - x1) - (y2 -y1)  def computeBCoeff(x1, x2, y1, y2, k):     """     >>> print(computeBCoeff(10., 42., 128., -64., 6.))     -384.0     """     return -k*(x2 - x1) + (y2 -y1)  def findLowerBound(xvalues, x):     """     >>> findLowerBound([10., 30.], 9.)     -1     >>> findLowerBound([10., 30.], 10.)     0     >>> findLowerBound([10., 30.], 15.)     0     >>> findLowerBound([10., 30.], 30.)     0     >>> findLowerBound([10., 30.], 31.)     1     >>> findLowerBound([10., 30., 40.], 9.)     -1     >>> findLowerBound([10., 30., 40.], 10.)     0     >>> findLowerBound([10., 30., 40.], 15.)     0     >>> findLowerBound([10., 30., 40.], 30.)     0     >>> findLowerBound([10., 30., 40.], 40.)     1     >>> findLowerBound([10., 30., 40.], 41.)     2     """     if xvalues[-1] == x:         return len(xvalues) - 2      left = bisect_left(xvalues, x)      if left >= len(xvalues):         return len(xvalues) - 1      if (xvalues[left]==x):         return 0 if left == 0 else left - 1      return left - 1 ``` 

Postgresql: percentile_cont interpolation of values

I have a couple of questions regarding percentile_cont in Postgresql. I would like to know which algorithm it uses for interpolating values when using percentile_cont ? My second question is why is it that we have to use order by ? It just seems weird to me that it doesn’t order values by default when using percentile_cont? I fail to see why this choice has been made. Would someone please explain the reason behind such a choice?