Problem with DifferentialEquations and overwriting methods

General Usage
1

I’ve been trying to integrate some simple parameterized ODEs using the DifferentialEquations.jl package; my normal practice from Matlab has been to define a general function with three inputs (time, state, parameters) and then use an anonymous function to specify parameters before feeding the anonymous function into an integrator.

In Julia (v0.5) I’ve tried the following and even though the function f seems to be redefined correctly, the integration results do not change.

using DifferentialEquations

du_dt(t, u, p) = p[1]*u[1]

f(t, u) = du_dt(t, u, [-1.])
assert(f(0, 1.) == -1.)
prob1 = ODEProblem(f, [1.], (0., 1.))
sol1 = solve(prob1)

f(t, u) = du_dt(t, u, [-2.])
assert(f(0, 1.) == -2.)
prob2 = ODEProblem(f, [1.], (0., 1.))
sol2 = solve(prob2)

sol1(1.)
sol2(1.)

I end up with sol1(1.) == sol2(1.) when they should be different since I’ve used a different f.

Am I misunderstanding something here or is there a bug? Is there a better pattern for doing this sort of thing?

Thanks
David

2

There are a few things here. First of all, the issue you’re having is that

f(t, u) = du_dt(t, u, [-1.])

is not an anonymous function, and when you overwrite a function like this, bad things can happen. You should instead use an anonymous function:

f = (t, u) -> du_dt(t, u, [-1.])

and you should be good. However, I’d like to lay out the whole story here to help you out a little bit more.

In-Place Functions

You might want to take another look at the bottom of the ODE tutorial. In DifferentialEquations.jl, the 3-parameter function is used for the “in-place” updating function. For example,

function lorenz(t,u,du)
 du[1] = 10.0(u[2]-u[1])
 du[2] = u[1]*(28.0-u[3]) - u[2]
 du[3] = u[1]*u[2] - (8/3)*u[3]
end

Because this mutates du as output instead of outputing vectors, it avoids allocations and is fast. The downside is that this will trip up some ex-MATLAB users who will write the “parameter function” as you did.

There are a few things you can do.

Closures

The simplest way of handling parameters is to use closures. What your “true” equation should be, in the in-place updating form, is:

du_dt(t, u, du, p) = (du[1] = p[1]*u[1])

To make this something the solver can handle, you use a closure which saves the parameter vector into the function:

f = (t,u,du) ->  du_dt(t, u, du, [2.0])

Now f(t,u,du) has the parameter information in it, and you can give it to the solver.

ParameterizedFunctions.jl

A much more interesting way to define parameters is via ParameterizedFunctions.jl. Using the provided macro, we can define the Lorenz equation using:

using ParameterizedFunctions
f = @ode_def_nohes Lorenz begin
  dx = σ*(y-x)
  dy = x*(ρ-z) - y
  dz = x*y - β*z
end σ=>10.0 ρ=>28.0 β=(8/3)

The macro automatically will transform this syntax into an in-place function, and use SymEngine to do symbolic calculations to both simplify the function and automatically calculate things like Jacobians and Inverse Jacobians (which the solvers will use to get better performance). So if your ODE is a system of ODEs, this is probably the preferred way to do this because not only is it pretty, but it also allows the solvers to hyper-optimize (it’s not entirely used right now, but that’s only because the dev API for that was just finished last night).

Note that if you defining the function using ParameterizedFunctions.jl is required for parameter sensitivity analysis, parameter estimation, and other functionality which requires explicit handling of parameters.

Manual ParamaterizedFunction

The macro shown above cannot handle every case. It’s great for systems of ODEs/SDEs, but for a general large system of ODEs (usually derived from a PDE) it won’t work. In that case, you’ll have to define the ParameterizedFunction manually. To do that, you use Julia’s call-overloading syntax. Here’s an example:

type  LotkaVolterra <: ParameterizedFunction
         a::Float64
         b::Float64
end
f = LotkaVolterra(0.0,0.0)
(p::LotkaVolterra)(t,u,du) = begin
         du[1] = p.a * u[1] - p.b * u[1]*u[2]
         du[2] = -3 * u[2] + u[1]*u[2]
end

Note that what you do is you make a type, and you overload the call on the type so that way “it looks like a function”, but it’s actually a type holding your parameters, and those parameters can then be used in the call. This syntax can be used for any system and will still work in parameter estimation routines and such.

Final Statement

While this is all still working, probably what’s under the most development is the documentation for the users. I hope this sets you (and others who find this) on the right path. Please let us know any feedback you have on documentation and usage.

2 Likes
3

Thanks, that’s very helpful. I realised that my function wasn’t anonymous but I didn’t appreciate the nasty things that can happen like this (though I did know about the dependant function issue).

I like the idea of the manual ParameterizedFunction (the simpler syntax one isn’t appropriate for my particular use case).

Thanks again
David

4

Dear Chris and Dawbarton, your post are very helpful, thank you. I use Python 3+ on daily basis, but would like to transfer my code to Julia. Very recently I tried to implement a short code for one of my student in Julia dealing with the fiiting of experimental data based on simple differential equation, namely dc/dt = -k1*c^k2. Please have look on the following code.

How can I start with parametrized ODE´s and data fitting in Julia?! Any suggestions will be very appreciated :slight_smile:

I had great troubles installing (still not working) DifferentialEquations on Windows 10 64-bit (JuliaPro 6.0.1 with mingw -32)

Thank you,
Pavel Dytrych

https://github.com/buchs-kevin/Pandas-Excel-Example/blob/master/Pandas-Excel-Example.py - very interesting

import numpy as np
import pandas as pd
from scipy import integrate, interpolate
from scipy import optimize
import matplotlib.pylab as plt

import matplotlib
matplotlib.rcParams[‘text.usetex’] = True
matplotlib.rcParams[‘text.latex.unicode’] = True

“”" Old - data given by hand
x_data = np.array([
0,
60,
120,
240,
360,
600,
1080,
2040,
3960,
7800
])

y_data = np.array([
970,
739,
699,
654,
609,
556,
482,
393,
254,
52
])
“”"

file = ‘data_kinetics.xlsx’
df = pd.read_excel(file, sheetname=‘List1’)
print(df)
print(“The list of row indicies”)
print(df.index)
print(“The column headings”)
print(df.columns)

x_data = df[‘t’]
y_data = df[‘y’]

x, y = x_data, y_data

def f_pseudo(y, t, k):
dy0 = -k[0]*y[0]**k[1]
return dy0

def f_pseudo2(y, t, k):
dy0 = -k[0]*(y[0]-k[1])**2
return dy0

def my_ls_func(x,teta):
“”“definition of function for LS fit
x gives evaluation points,
teta is an array of parameters to be varied for fit”“”
# create an alias to f which passes the optional params
f2 = lambda y,t: f_pseudo(y, t, teta)
# calculate ode solution, retuen values for each entry of “x”
r = integrate.odeint(f2,y0,x)
#in this case, we only need one of the dependent variable values
return r[:,0]

def f_resid(p):
“”" function to pass to optimize.leastsq
The routine will square and sum the values returned by
this function"“”
return y_data-my_ls_func(x_data,p)

#solve the system - the solution is in variable c
guess = [0,1] #initial guess for params
y0 = [970] #inital conditions for ODEs
(c,kvg) = optimize.leastsq(f_resid,guess) #get params
popt, pcov = optimize.leastsq(f_resid,guess)

Parameters for pseudo 1. order

print (c[0], c[1])
pf11=c[0]
pf12=c[1]

def my_ls_func2(x,teta):
“”“definition of function for LS fit
x gives evaluation points,
teta is an array of parameters to be varied for fit”“”
# create an alias to f which passes the optional params
f3 = lambda y,t: f_pseudo2(y, t, teta)
# calculate ode solution, retuen values for each entry of “x”
r = integrate.odeint(f3,y02,x)
#in this case, we only need one of the dependent variable values
return r[:,0]

def f_resid2(p):
“”" function to pass to optimize.leastsq
The routine will square and sum the values returned by
this function"“”
return y_data-my_ls_func2(x_data,p)

#solve the system - the solution is in variable c
guess2 = [0,1] #initial guess for params
y02 = [970] #inital conditions for ODEs
(c2,kvg) = optimize.leastsq(f_resid2,guess2) #get params
popt2, pcov2 = optimize.leastsq(f_resid2,guess2)

Parameters for pseudo 2. order

print (c2[0], c2[1])
ps21=c2[0]
ps22=c2[1]

xeval=np.linspace(min(x_data), max(x_data),500)
gls = interpolate.UnivariateSpline(xeval, my_ls_func(xeval,c), k=5, s=0)
gls2 = interpolate.UnivariateSpline(xeval, my_ls_func2(xeval,c2), k=5, s=0)

f_experimental = open(“Data_experimental.csv”, “w”)
for i in range(len(x_data)):
f_experimental.write(“{}, {} \n”.format(x_data[i], y_data[i]))
f_experimental.close()

f_data = open(“Data_fitting.csv”, “w”)
for i in range(len(xeval)):
f_data.write(“{}, {}, {} \n”.format(xeval[i], my_ls_func(xeval,c)[i], my_ls_func2(xeval,c2)[i]))
f_data.close()

f_L = open(“Fitting_parameters_pseudo_first_order.csv”, “w”)
f_L.write(“{}, {} \n”.format(pf11, pf12))
f_L.close()

f_F = open(“Fitting_parameters_pseudo_second_order.csv”, “w”)
f_F.write(“{}, {} \n”.format(ps21, ps22))
f_F.close()

plt.figure(1, figsize=(8, 6))
plt.xlabel(r’t’, fontsize=18)
plt.ylabel(r’q_t’, fontsize=18)
plt.title(r"Experimental data adsorption",fontsize=18)
plt.xticks(fontsize=18)
plt.yticks(fontsize=18)

plt.plot(x_data,y_data,‘*’, label=‘Experimental data’)
plt.plot(xeval, gls(xeval), ‘b-’, label=‘Pseudo-first order’)
plt.plot(xeval, gls2(xeval), ‘r-’, label=‘Pseudo-second order’)
legend = plt.legend(loc=‘best’, shadow=True, fontsize=‘x-large’)
legend.get_frame().set_facecolor(‘#FFFFFF’)
plt.savefig(“Rozklad_silanu.png”, dpi=300)
plt.show()

5

Hi @Pavel_Dytrych, could you please open up a different thread? You probably shouldn’t be hijacking a thread like this.

For parameter estimation, check out the documentation page:

http://docs.juliadiffeq.org/latest/analysis/parameter_estimation.html

As for JuliaPro, don’t use JuliaPro if you plan on adding packages. There are a lot of Julia packages, DifferentialEquations.jl included, which are simply not compatible with JuliaPro. I recommend just using a standard Julia installation instead, at least until this is fixed.

2 Likes