CODEWITHSUNDEEP
pythonscipydifferential-equations
Andrej
Andrej

Reputation: 437

How to properly implement scipy.integrate.Radau?

(I'll apreciate even a link to an example, so far I haven't found any.)

I am trying to use Radau from scipy.integrate to solve second order differential equation. For now I am trying just a simple example, so that I can uderstand how it works (unsuccessfully so far).

Let's say that I have following equation:

d^2y/dt^2 = C,

which means that y = (Ct^2)/2 + Bt.

Let's say, for example, y(1) = 1 and C = 2. Let's say that I want to find value of y for t = 10.

This is my code:

from scipy.integrate import Radau
import numpy as np

C = 2.0
y = np.zeros(2)

def fun(t, y):
  #y[0] = C*t
  y[1] = C
  return y

t0 = 1.0
y0 = np.array([1., 1.])
t_bound = 10.0

eq = Radau(fun, t0, y0, t_bound)
print(eq.n)

while(True):
  print(eq.t)
  print(eq.y)
  print(eq.status)
  if eq.status == 'finished':
    break
  eq.step()

Outputs is wrong. (If I uncomment that one commented line in fun definition, it also gives wrong answer. But I think that I shouldn't even have to tell solver that, right? I don't know this value usually.)

I think my biggest problem is that I am not really sure what should be passed as fun. Documentation says it should be right-hand side of the system, so I thought that first derivation should be in y[0], second derivation in y[1] etc.

What am I doing wrong? How should this be implemented?

Upvotes: 1

Views: 1201

Answers (1)

Lutz Lehmann
Lutz Lehmann

Reputation: 25992

Atm you are solving 

y0' = y0,  y0(1)=1
y1' = 2,   y1(1)=1

which has the solution y0(t)=exp(t-1) and y1(t)=2*t-1 which is certainly not what you want. You want the first order system

y0' = y1
y1' = C

so that you need

def fun(t,y): return [y[1], C]

Then the solution is y1(t)=C*t+B=2*t-1 and y0(t)=0.5*C*t^2+B*t+A=t^2-t+1 and the integration ends correctly with eq.y = [91. 19.].

Upvotes: 3

Related Questions