External polynomial supply to DifferentialEquation.jl

General Usage
1

Hi everyone

I have a problem solving the Riccati type differential equation.

\frac{dU}{dt} = \delta \frac{dD}{dt} - \lambda U - \gamma U^2

I have an array of dimension M for D from which, in my opinion we can get a polynomial and then supply the polynomial into DifferentialEquations.jl routine inside which we can then differentiate and convert back to an array of dimension M.

I also tried using the finite difference approach, but the right-hand side with \frac{dD}{dt} has dimension M-1, which is a problem.

Any help is appreciated.

2

Use DataInterpolations.jl and grab its derivative?

3

Thank you so much for your quick response.

I have used the DataInterpolations.jl and generated a p=CubicSpline(D,t) but taking the derivative with the command derivative(p) is not working. To make it work, I have to use convert(Polynomial,p), and then the differentiation is possible, but the nodes are of magnitude 1e16.

Experiment:

  1. Considering M=11 nodes, with t=collect(0.0:1.0:10) and u=ones(M)
  2. p=CubicSpline(u,t)
  3. c=convert(Polynomial,p)
  4. p_2 = derivative(c)
  5. p_2.(t) = [1.0, 1045.0, 3.81707e8,1.0225e12 ,2.93529e14, 2.41548e16 ,8.95379e17,1.90888e19,2.71108e20,2.82139e21,2.29657e22]

Apart from this being of large magnitude, it is working.

Is there any possibility that we can have a polynomial of a certain order (like maybe 5 or something) but with a number of nodes that is fixed like M in this case?

4

The problem is the conversation to a Polynomial. You could use this instead:

using DataInterpolations

M = 11
ts = 0.0:10.0
u = ones(M)

p = CubicSpline(u, ts)
D = [ DataInterpolations.derivative(p, t) for t in ts]

If you want the derivative as a function, you could use ForwardDiff or simply define

dp_dt = t -> DataInterpolations.derivative(p, t)

and use that in your ODE problem. (But make sure that p in this setting is not a global variable, but a parameter of the ODE, otherwise it will be unnecessarily slow.)

PS: By the way, you can embed source code with one or three backticks, i.e. ``` source code ```.

5

Thanks, @SteffenPL, and @ChrisRackauckas. I was able to make it work through your instructions. I was wondering about what if we use Chebyshev polynomial instead of CubicSpline. Would it work provided this number of nodes would go up and polynomial would of order M?

using Polynomials 
M = 11
u = ones(M)
p = ChebyshevT(u)
p1 = convert(Polynomial, c)

I would like to have a Chebyshev polynomial of order, say 5, instead of going through all M nodes.