How to solve a differential equation using Python's built-in odeint function?

Question

I want to solve these differential equations with given initial conditions:

(3x-1)y''-(3x+2)y'+(6x-8)y=0, y(0)=2, y'(0)=3 

ans should be

y=2*exp(2*x)-x*exp(-x)

here is my code:

 def g(y,x): y0 = y[0] y1 = y[1] y2 = (6*x-8)*y0/(3*x-1)+(3*x+2)*y1/(3*x-1) return [y1,y2] init = [2.0, 3.0] x=np.linspace(-2,2,100) sol=spi.odeint(g,init,x) plt.plot(x,sol[:,0]) plt.show() 

but what I get is different from the answer. What did I do wrong?

Answer 1

Here are a few things wrong. First, your equation seems to be

(3x-1) y '' - (3x + 2) y '- (6x-8) y = 0; y (0) = 2, y '(0) = 3

(note the sign of the term in y). For this equation, your analytical solution and definition of y2 are correct.

Secondly, as @Warren Weckesser says, you must pass 2 parameters as y to g : y[0] (y), y[1] (y ') and return your derivatives y' and y '.

Thirdly, your initial conditions are set for x = 0, but your x-grid for integration starts at -2. From the docs for odeint this parameter t in their call signature description:

odeint(func, y0, t, args=(),...) :

t: array A sequence of time points to solve for y. The initial value must be the first element of this sequence.

So, you must integrate the beginning at 0 or provide initial conditions starting from -2.

Finally, your integration range covers the feature at x = 1/3. odeint may have a bad time here (but apparently not).

Here's one approach that seems to work:

 import numpy as np import scipy as sp from scipy.integrate import odeint import matplotlib.pyplot as plt def g(y, x): y0 = y[0] y1 = y[1] y2 = ((3*x+2)*y1 + (6*x-8)*y0)/(3*x-1) return y1, y2 # Initial conditions on y, y' at x=0 init = 2.0, 3.0 # First integrate from 0 to 2 x = np.linspace(0,2,100) sol=odeint(g, init, x) # Then integrate from 0 to -2 plt.plot(x, sol[:,0], color='b') x = np.linspace(0,-2,100) sol=odeint(g, init, x) plt.plot(x, sol[:,0], color='b') # The analytical answer in red dots exact_x = np.linspace(-2,2,10) exact_y = 2*np.exp(2*exact_x)-exact_x*np.exp(-exact_x) plt.plot(exact_x,exact_y, 'o', color='r', label='exact') plt.legend() plt.show()