Help with Octave for system of ODEs


by McLaren Rulez
Tags: octave, odes
McLaren Rulez
McLaren Rulez is offline
#1
Dec10-12, 02:48 AM
P: 261
Hi,

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]
[tex]\frac{db}{dt}=1[/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);
xdot=zeros(2,1)
xdot(1)=2*x(1)*x(2)
xdot(2)=1
endfunction

t=linspace(0,10,100);
x=lsode("f",[1;0],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.
Phys.Org News Partner Science news on Phys.org
SensaBubble: It's a bubble, but not as we know it (w/ video)
The hemihelix: Scientists discover a new shape using rubber bands (w/ video)
Microbes provide insights into evolution of human language
gsal
gsal is offline
#2
Dec10-12, 10:23 PM
P: 838
Well...don't know ODEs nor Octave, but a quick look at python's scipy revealed a similar function.

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)
ax1.plot(t,x[:,0])
ax1.set_ylabel('a(t)')

ax2 = fig.add_subplot(212)
ax2.plot(t,x[:,1])
ax2.set_xlabel('t')
ax2.set_ylabel('b(t)')

plt.show()
See attached plot, too.
Attached Thumbnails
ab-ode.png  


Register to reply

Related Discussions
Transforming a system of PDEs into a first order system of ODEs Calculus & Beyond Homework 0
System of ODEs Calculus & Beyond Homework 5
A System of Three ODEs Calculus & Beyond Homework 1
From a System of 1st ODE to a 2nd ODE and back to the system of 1st ODEs Differential Equations 1
system of ODEs Differential Equations 2