1

I'm having a little trouble trying to understand what's wrong with me code, any help would be extremely helpful. I wanted to solve this simple equation

enter image description here

However, the values my code gives doesn't match with my book ones or wolfram ones as y goes up as x grows.

import matplotlib.pyplot as plt
from numpy import exp
from scipy.integrate import ode

# initial values 
y0, t0 = [1.0], 0.0

def f(t, y):
    f = [3.0*y[0] - 4.0/exp(t)]
    return f

# initialize the 4th order Runge-Kutta solver
r = ode(f).set_integrator('dopri5')
r.set_initial_value(y0, t0)

t1 = 10
dt = 0.1
x, y = [], []
while r.successful() and r.t < t1:
    x.append(r.t+dt); y.append(r.integrate(r.t+dt))
    print(r.t+dt, r.integrate(r.t+dt))
Orbitz
  • 53
  • 5
  • 1
    Try a smaller `dt`. You'll notice you get farther before it blows up. When you are required to take excessively small steps when the exact solution is quite smooth, you should suspect that your ODE is "stiff": https://en.wikipedia.org/wiki/Stiff_equation – Klimaat Apr 18 '17 at 21:25
  • Oh, I see! With smaller dt the equation works fine, now I understood, I thought I was doing something wrong because it's the first time I use this ode function. Well, thanks. – Orbitz Apr 18 '17 at 22:02
  • @klimaat. Could you please post your comment as an answer so OP can accept it and the question can be marked as answered? Also, you will both get a little reputation out of it which is nice. – Mad Physicist Apr 19 '17 at 04:31

1 Answers1

1

Your equation in general has the solution

y(x) = (y0-1)*exp(3*x) + exp(-x)

Due to the choice of initial conditions, the exact solution does not contain the growing component of the first term. However, small perturbations due to discretization and floating point errors will generate a non-zero coefficient in the growing term. Now at the end of the integration interval this random coefficient is multiplied by exp(3*10)=1.107e+13 which will magnify small discretization errors of size 1e-7 to contributions in the result of size 1e+6 as observed when running the original code.

You can force the integrator to be more precise in its internal steps without reducing the output step size dt by setting error thresholds like in

r = ode(f).set_integrator('dopri5', atol=1e-16, rtol=1e-20)

However, you can not avoid the deterioration of the result completely as the floating point errors of size 1e-16 get magnified to global error contributions of size 1e-3.

Also, you should notice that each call of r.integrate(r.t+dt) will advance the integrator by dt so that the stored array and the printed values are in lock-step. If you want to just print the current state of the integrator use 

    print(r.t,r.y,yexact(r.t,y0))

where the last is to compare to the exact solution which is, as already said,

def yexact(x,y0): 
    return [ (y0[0]-1)*exp(3*x)+exp(-x) ]
Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51