Register to reply

Help with Octave for system of ODEs

by McLaren Rulez
Tags: octave, odes
Share this thread:
McLaren Rulez
#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
Scientists develop 'electronic nose' for rapid detection of C. diff infection
Why plants in the office make us more productive
Tesla Motors dealing as states play factory poker
gsal
#2
Dec10-12, 10:23 PM
P: 890
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