Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

Help with Octave for system of ODEs

  1. Dec 10, 2012 #1

    I am having a lot of trouble with Octave as I try to solve a system of ODEs. Any help is appreciated, I am a complete newbie with Octave and numerical solving.

    Let's try a very simple one. Suppose I had a pair of ODEs with a and b being functions of time

    [tex] \frac{da}{dt}=2ba[/tex]
    Initial conditions are a(0)=1, b(0)=0

    This is clearly the solved by [tex]a(t)=e^{t^{2}}[/tex] [tex]b(t)=t[/tex] My Octave code was this:

    function xdot=f(x,t);


    I want to plot a(t) against t or b(t) or some combination of a and b against t. Here are my issues

    1) The t=linspace() part. What numbers are appropriate? Sometimes, I got an error saying convergence failure but this combinations worked through blind luck. In general, what should I choose and why does it seem to matter? As I understand, this tells Octave to take t from 0 to 10 and have 100 intervals. I thought any numbers there would have worked?

    2) This is more important. I tried plot(t,x(1)) but I got a blank plot. plot(t,x(2)) also gave me a blank plot. plot(t,x) gave me something but it's really weird. Isn't x now a column vector? I'm not sure what exactly lsode outputs here. What should be the correct command to get a(t) against t, which must of course be an exponential t squared against t graph?

    There's also the fact that when I do it for my actual set of ODEs which are slightly more complicated, it inevitably hits an error or gets something 'x' undefined at a certain column and certain line. I'm quite lost :(

    Thank you for you help.
    Last edited: Dec 10, 2012
  2. jcsd
  3. Dec 10, 2012 #2
    Well...don't know ODEs nor Octave, but a quick look at python's scipy revealed a similar function.

    Code (Text):

    import numpy as np
    from scipy.integrate import odeint

    import matplotlib.pyplot as plt

    def dxdt(x,t):
        xdot = np.zeros(2)
        xdot[0] = 2.0*x[0]*x[1]
        xdot[1] = 1.0
        return xdot
    t = np.linspace(0,5,51)
    x = odeint(dxdt, [1.0,0.0], t)

    fig = plt.figure()

    ax1 = fig.add_subplot(211)

    ax2 = fig.add_subplot(212)


    See attached plot, too.

    Attached Files:

Share this great discussion with others via Reddit, Google+, Twitter, or Facebook