TransWikia.com

ODE coupled to timedependent polynomial

Mathematica Asked by Running_Quark on June 6, 2021

I’m trying to solve a differential equation a'[t]== f[a[t], y[t]], where for every time step y[t] is a real solution to a fourth-order polynomial p[a[t],y[t]]. As an old Python programmer, it is trivial to do with a loop and a step method like Runge-Kutta 4. But is there a nice, perhaps quicker, way to do this with DSolve/NDSolve or similar? I hope there is a fast and simple way to do it in Mathematica.

As an example:

ODE:

-0.7-0.3/a[t]^3-50/3 (-1+y[t]) (2+(2-y[t]) y[t])+a'[t]^2/a[t]^2==0

Polynomial:

(-82.6333+0.3/a[t]^3) y[t]+50 y[t]^2+49.3 y[t]^3-(50 y[t]^4)/3==0

The initial conditions are not super-relevant and a[1]==1 is sufficient.

2 Answers

First we have to decide which initial conditions to use for the DAE system.

soly=Values@Solve[ (-82.6333 + 0.3/a[t]^3) y[t] + 50 y[t]^2 + 49.3 y[t]^3 - (50 y[t]^4)/3 == 0, y[t]] /. a[t] -> 1 // Chop
(*{{0}, {-1.45186}, {3.41291}, {0.996957}}*)

OP looks for a branch 0<y<1 (see comment)

y1=soly[[-1]][[1]] (*0.996957*)

Now we solve the DAE-system with intial conditions a[1]==1, y[1]==y1

sola = NDSolveValue[{-0.7 - 0.3/a[t]^3 -50/3 (-1 + y[t]) (2 + (2 - y[t]) y[t]) + a'[t]^2/a[t]^2 ==0
, (-82.6333 + 0.3/a[t]^3) y[t] + 50 y[t]^2 +49.3 y[t]^3 - (50 y[t]^4)/3 == 0 , a[1] == 1, y[1] == y1},a, {t, 1, 2}]
Plot[sola[t] , {t, 1, 2}]

enter image description here

Answered by Ulrich Neumann on June 6, 2021

Note: I am not very experienced with NDSolve and differential algebraic equations, so the following solution is probably not the nicest one.


Normally, NDSolve is able to handle differential algebraic equations automatically. However, such equations usually have more than one solution, and we need to do some manual work to select the one we need. While the below works, I am also interested in the best way to do this.

Taking your equations,

deq = -0.7 - 0.3/a[t]^3 - 50/3 (-1 + y[t]) (2 + (2 - y[t]) y[t]) + 
   a'[t]^2/a[t]^2 == 0

eq = (-82.6333 + 0.3/a[t]^3) y[t] + 50 y[t]^2 + 
   49.3 y[t]^3 - (50 y[t]^4)/3 == 0

we can try a direct solution:

sol = NDSolve[
  {deq, eq, a[0] == 1},
  {a, y},
  {t, 0, 10}
  ]

But there is a warning that the computed initial conditions do not match the initial conditions we provided manually:

enter image description here

Let's verify them:

In[4]:= {a[0], y[0]} /. sol
Out[4]= {{2.7121*10^6, 1.}}

Indeed, a[0] is far from 1. I assume (but I don't know) that this might be due to Mathematica choosing the "wrong" y[0] solution to eq.

Let us do some manual investigation then. eq is a fourth-order polynomial, so for a given a[0] we have four possible y[0] values. But not all real y[0] values are acceptable, as some of them lead to a negative a'[0] when solving deq.

This plot shows the range of acceptable y[0] values:

Plot[
 ad /. Solve[deq /. t -> 0 /. a[0] -> 1 /. a'[0]^2 -> ad, ad] // 
  Evaluate,
 {y[0], -3, 4}
 ]

enter image description here

For a[0] == 1, we get these possible y[0] (all real):

ysol = Solve[eq /. t -> 0 /. a[0] -> 1]
(* {{y[0] -> -1.45186}, {y[0] -> 0.}, {y[0] -> 0.996957}, {y[0] -> 3.41291}} *)

But only the first and third one give valid a'[0] values:

NSolve /@ (deq /. t -> 0 /. a[0] -> 1 /. ysol)
(* {{{Derivative[1][a][0] -> -11.1386}, {Derivative[1][a][0] -> 
    11.1386}}, {{Derivative[1][a][0] -> 
    0. - 5.68624 I}, {Derivative[1][a][0] -> 
    0. + 5.68624 I}}, {{Derivative[1][a][
     0] -> -0.920791}, {Derivative[1][a][0] -> 
    0.920791}}, {{Derivative[1][a][0] -> 
    0. - 10.6062 I}, {Derivative[1][a][0] -> 0. + 10.6062 I}}} *)

enter image description here

We can compute a solution e.g. using the third one:

NDSolve[
 {deq, eq, a[0] == 1, ysol[[3]] /. Rule -> Equal},
 {a, y},
 {t, 0, 10}
 ]

(EDIT: I originally missed OP's comment that y[0] should be in [0,1]. It seems that it is indeed the 3rd solution that satisfies this.)

Note that this is still not the only solution for the given y[0], as a'[0] can be positive or negative. You can compute a desired a'[0] and add that to the initial conditions as well.


Now my personal question to other people who might write an answer:

Is there a simpler way to choose our desired solution? Is there a way to avoid manually computing the initial conditions and plugging them in? For example, can we somehow just specify that we want a'[0] >= 0 instead of computing an explicit value for a'[0] by hand and plugging that in? Can we do the same asking that 0 < y[0] < 1?

Answered by Szabolcs on June 6, 2021

Add your own answers!

Ask a Question

Get help from others!

© 2024 TransWikia.com. All rights reserved. Sites we Love: PCI Database, UKBizDB, Menu Kuliner, Sharing RPP