2

i am a newbie to python. I have a simple differential systems, which consists of two variables and two differential equations and initial conditions x0=1, y0=2:

dx/dt=6*y
dy/dt=(2t-3x)/4y

now i am trying to solve these two differential equations and i choose odeint. Here is my code:

import matplotlib.pyplot as pl
import numpy as np
from scipy.integrate import odeint

def func(z,b):
    x, y=z
    return [6*y, (b-3*x)/(4*y)]    

z0=[1,2]
t = np.linspace(0,10,11)
b=2*t
xx=odeint(func, z0, b)
pl.figure(1)
pl.plot(t, xx[:,0])
pl.legend()
pl.show()

but the result is incorrect and there is a error message:

enter image description here

Excess work done on this call (perhaps wrong Dfun type).
Run with full_output = 1 to get quantitative information.

I don't know what is wrong with my code and how can i solve it. Any help will be a useful to me.

percusse
  • 2,808
  • 1
  • 12
  • 26
HH624
  • 35
  • 1
  • 5
  • 1
    Same error message and context as in http://stackoverflow.com/questions/34616775/using-odeint-to-solve-two-first-order-differential-equations-with-a-vektor-as-in. Just a simpler function. – hpaulj Jan 05 '16 at 19:52
  • 1
    Possible answer in http://stackoverflow.com/questions/28092456/scipy-odeint-giving-lsoda-warning. This talks about the same `lsoda` warning that I get running your function. He's saying the problem isn't with the code, it's with the behavior of the ode. – hpaulj Jan 05 '16 at 20:06
  • thank u, i am reading the question of the links. i want to know, have i used the odeint right to solve these two differential equations? i have never used odeint for solving two differential equations – HH624 Jan 05 '16 at 20:18
  • You used the interface right, however you did not implement the announced ODE. There is a difference in solving `y'(t)=f(2t,y(t))` and `y'(t)=f(t,y(t))` that is not resolved by doubling the sample points. – Lutz Lehmann Jan 05 '16 at 20:27
  • `integrate.odeint(func, [1.,2.],np.linspace(0,.6,10))` runs fine. But when I up the range to .7 I start getting this warning. As `LutzL` notes, `y` is getting too close to 0, and your derivative is blowing up. – hpaulj Jan 05 '16 at 22:40

1 Answers1

5

Apply trick to desingularize the division by y, print all ODE function evaluations, plot both components, and use the right differential equation with the modified code

import matplotlib.pyplot as pl
import numpy as np
from scipy.integrate import odeint

def func(z,t):
    x, y=z
    print t,z
    return [6*y, (2*t-3*x)*y/(4*y**2+1e-12)]    

z0=[1,2]
t = np.linspace(0,1,501)
xx=odeint(func, z0, t)
pl.figure(1)
pl.plot(t, xx[:,0],t,xx[:,1])
pl.legend()
pl.show()

and you see that at t=0.64230232515 the singularity of y=0 is assumed, where y behaves like a square root function at its apex. There is no way to cross that singularity, as the slope of y goes to infinity. At this point, the solution is no longer continuously differentiable, and thus this is the extremal point of the solution. The constant continuation is an artifact of the desingularization, not a valid solution.

graph of x and y over t

Lutz Lehmann
  • 21,318
  • 2
  • 18
  • 42