Help on defining function for ODE solver

General Usage
1

Newbie here. I am trying to write ODE function with vector calculation in it. And I started with very basic examples. I tried one with expression definition and another with function definition. What puzzles me is that they show different outputs. The function definition won’t show any change of output. I am attaching the code and results below.

Thanks in advance!
Wenzhe

using DifferentialEquations

function g!(dx,x,p,t)
A = p
dx = A*x
return dx
end

A=rand(2,2)
u=rand(2,1)
g(x,p,t)=A*x

prob=ODEProblem(g,u,(0.0,10.0),A)
sol=solve(prob)
println(sol)

prob=ODEProblem(g!,u,(0.0,10.0),A)
sol1=solve(prob)
println(sol1)

Interpolation: Automatic order switching interpolation
t: [0.0, 0.08492876287055985, 0.2632053877829872, 0.4919846968605104, 0.7721098086852171, 1.1031553910836651, 1.4829484225702116, 1.9084842676628948, 2.3761019883528935, 2.8818659161037585, 3.421794538745908, 3.9920109522139495, 4.588871868243653, 5.209011997880401, 5.849382858303235, 6.507234484776225, 7.180174894155538, 7.866060691297262, 8.563041192451601, 9.269505209830445, 9.984062996135261, 10.0]
u: [[0.8361627424041629; 0.993017552066968], [0.9472638967939113; 1.1317511689403286], [1.2321309385945562; 1.4877445910570049], [1.7294965772234105; 2.1098827839335264], [2.624588873507298; 3.2304929303957604], [4.304689172194103; 5.335381705725871], [7.606080961094069; 9.473643810771797], [14.41076700026382; 18.006252708395813], [29.108715532142092; 36.44048598309004], [62.30095382851463; 78.07554925349324], [140.4172088030801; 176.0681524042761], [331.2930799027888; 415.51972984293917], [813.6860331122462; 1020.686041610875], [2069.8356283051935; 2596.548930391782], [5428.089957953568; 6809.555143923475], [14614.634122998978; 18334.299617345787], [40252.47734924911; 50497.620521076715], [113046.93821852087; 141820.12961063726], [322831.7508407441; 405000.5566720919], [935172.0435733753; 1.1731971853211848e6], [2.7421810613658247e6; 3.4401364540926614e6], [2.808777742759984e6; 3.523683703285686e6]]
retcode: Success
Interpolation: Automatic order switching interpolation
t: [0.0, 1.0e-6, 1.1e-5, 0.00011099999999999999, 0.0011109999999999998, 0.011110999999999996, 0.11111099999999996, 1.1111109999999995, 10.0]
u: [[0.8361627424041629; 0.993017552066968], [0.8361627424041629; 0.993017552066968], [0.8361627424041629; 0.993017552066968], [0.8361627424041629; 0.993017552066968], [0.8361627424041629; 0.993017552066968], [0.8361627424041629; 0.993017552066968], [0.8361627424041629; 0.993017552066968], [0.8361627424041629; 0.993017552066968], [0.8361627424041629; 0.993017552066968]]

2

This one is supposed to be a mutating function, so instead of creating a new vector it needs to modify dx, i.e.:

function g!(dx,x,p,t)
  A = p
  mul!(dx,A,x)
  return nothing
end

or the easier (but slower) way:

function g!(dx,x,p,t)
  A = p
  dx .= A*x
  return nothing
end

Notice how it can work when returning nothing: the output is ignored in the mutating form.

2 Likes
3

Your suggestion works! Thank you very much!

The key thing seems to be changing dx = Ax to dx.=Ax. Could you educate me a little bit about the difference between the two? Does the first one created a new variable?

Thanks

4

This tutorial probably will help in multiple ways: Optimizing DiffEq Code . It describes the in-place usage and how to write non-allocating code, and demonstrates mutation and .= as part of that.

2 Likes