Having fun with Ordinary Differential Equations in Clojure using com.konato.ode part 2

Introduction

In the first tutorial, we introduced the com.konato.ode library to simulate differential equations with a classical damped spring mechanical system. We used as a source the system of differential equations. In this second tutorial, we will show how to use the functions transfer2ode and state2ode to convert the same system presented on the form of transfert function and state space into systems of differential equations and the simulate them using the odeslolve function.

Converting Transfert functions

Let say, we want to convert a simple integrator with a fixed input of 1. The transfer function of an integrator is 1/s. The function transfer2ode require three paramaters, the numerator, the denominator and the input function. The numerator should be a least one degree less than the denominator.

user> (def fps (transfer2ode [1] [1 0] [1]))
#'user/fps
user> fps
{:dotfps {:x1 #<ode$tf2odedotlast__2126$fn__2131 com.konato.ode$tf2odedotlast__2126$fn__2131@4d40df>}, :outfps {:y1 #<ode$tf2odeout__2157$fn__2159 com.konato.ode$tf2odeout__2157$fn__2159@1de041e>}}

fps is an hashmap who contains differential equation and the output equation. We can now use them to simulate our system for a period t of 2 seconds using odesolve function.

user> (keepdecim (last (:y1 (odesolve rk4 (:dotfps fps) {:x1 0.0} :t 0 2 0.01 (:outfps fps)))) 5)
2.0

This is the result of the integration of 1 over a period t 2 second as expected.

Now, we will use the transfert function the classical damped spring mechanical system used in part 1:

Transfer function of a classical damped spring oscillator

Transfer function of a classical damped spring oscillator

We will create a function who can simulate the system to  allow validation with the results obtained in part 1:

(defn transferdampedho [m fv k xini xdotini T0 TF DT]
"Test the Runge-Kutto order 4 method and the transfer2ode for a classic damped harmonic oscilator."
 (let [
 DEN [1 (/ fv m) (/ k m)]
 NUM [0 1]
 inp [0]
 fps (transfer2ode NUM DEN inp)
 ]
 (odesolve rk4 (:dotfps fps) {:x2 0.0, :x1 0.9} :t T0 TF DT (:outfps fps))
))

Then we can study the error between the classical versus the transfer function conversion:

user> (= 0 (reduce + (map #(- %1 %2) (:x (dampedho 1 1 1 1.0 0.0 0.0 2.0 0.01)) (:y1 (transferdampedho 1 1 1 1.0 0.0 0.0 2.0 0.01)))))
true

The result are the same as expected.

Converting State Space

If we have a State Space representation of a system, we can use the state2ode function to convert it into differentials and outputs equation who then can be used by odesolve. Let say, we have our classical damped spring oscillator in the form of a State Space:

dampedspringstatesmall

We will create a function who can simulate the system to  allow validation with the results obtained in part 1:

(defn statedampedho [m fv k xini xdotini T0 TF DT]
"Test the Runge-Kutto order 4 method and the state2ode for a classic damped harmonic oscilator."
 (let [
 A [[0 1][(* -1 (/ k m)) (* -1 (/ fv m))]]  
 B [[0][1]]
 C [[1 0]]
 D [[0]]
 U [0]
 fps (state2ode A B C D U)
 ]
 (odesolve rk4 (:dotfps fps) {:x2 0.0, :x1 0.9} :t T0 TF DT (:outfps fps))
))

Then we can study the error between the classical versus the transfer function conversion:

user>  (= 0 (reduce +(map #(- %1 %2) (:x (dampedho 1 1 1 1.0 0.0 0.0 2.0 0.01)) (:y1 (statedampedho 1 1 1 1.0 0.0 0.0 2.0 0.01)))))
true

The result are the same as expected.

Conclusion

In this second part, we studied usage of the conversion function included in the com.konato.ode library. Conversion of Transfer functions and State Space were demonstrated. Subsequent tutorials will show more possibilities of the library.

2 Comments

Having fun with Ordinary Differential Equations in Clojure using com.konato.ode part 1

Introduction

This tutorial show how to use the com.konato.ode library to simulate a damped spring system in Clojure. Clojure is a relatively new functional langage using a java virtual machine created by Rich Hickey. At Konato, we wanted to explore this language to see if we can use it as a development language for our products. We have chosed a non trivial project to be used to validate the use of Clojure. An ordinary differential equations solvers and simulations library looked like an excellent challenge.  We hope you will have fun using it either as an eductationnal tool or to simulate real systems.

A classical damped harmonic oscillator has been chosen to introduce the usage of the library. In following tutorials, it will be shown how to use State Space or Transfer function instead of the differential equations,  control laws will be demonstrated.

The system and his equation

We have a mass M attached to a spring of coefficient k and a coefficient of friction Fv. The initials conditions are mass is at position x=0.9 with a speed 0 and 0 force applied to the mass.

A classical damped spring with friction

The following differential equation describe the system:

diffequdampedspringsmall

Using the simulation library

First make sure the com_konato_ode.tar archive or the source of the library is in your java classpath.

To use the library, you  have to enter the equations using anonymous functions syntax in a Clojure map:

       dotfps {:x #(:xdot %), :xdot #(/ (+ (* -1.0 k (:x %)) (* -1.0 fv (:xdot %))) m)}
       outfps {:x #(:x %), :t #(:t %)}

dotfps are the dotted equations to solve, keywords represent the resolved variables after integration. A dotted equation is required for each order a the system. Then you select the desired outputs via outfps. Simulation are done using the odesolve functions where you specify the odesolver, the initials conditions, time variable, start time, final time, interval when using a fixed step solver and the output functions.

(odesolve rk4 fps {:xdot 0.0, :x 0.9} :t T0 TF DT outfps)

Here is a complete function who allows mass, spring coefficent k and damping fv to be specified for a simulation with initials conditions position x  = 0.9 and speed xdot = 0:

(use 'com.konato.ode)
(defn dampedho [m fv k xini xdotini T0 TF DT]
"Test the Runge-Kutta order 4 method for a classic damped harmonic oscilator."
 (let [
       fps {:x #(:xdot %), :xdot #(/ (+ (* -1.0 k (:x %)) (* -1.0 fv (:xdot %))) m)}
       outfps {:x #(:x %), :t #(:t %)}
      ]
      (odesolve rk4 fps {:xdot 0.0, :x 0.9} :t T0 TF DT outfps)
))

Simulations can be done now easily for various roots:

  (dampedho 1 25 2 1.0 0.0 0.0 2.0 0.01); Real and unequal roots
  (dampedho 1 1 1 1.0 0.0 0.0 2.0 0.01); Complex roots
  (dampedho 1 2 100 0.9 0.0 0.0 2.0 0.01); Complex roots from G.A. Kornm,1989, Interactive Dynamic System Simulation, example 1.4 
(dampedho 1 4 4 1.0 0.0 0.0 2.0 0.01); Real and equal roots

Exporting the data to a file

I have created a small function who allow to export on a file results of simulations in a row based output. External software like gnuplot can then be used to visualize the data produced.

To save the position x and the time t, evaluate:

(savedata "datafile1.dat" (:x res) (:t res))

Data is then saved in a file, ready to be processed and plotted using gnuplot.

Plotting with gnuplot

If you are on Linux and have gnuplot installed, you can adapt the following script to plot your results:

#!/bin/sh
gnuplot << EOF
set terminal png large enhance
set output "exdampfv2k100.png"
set xlabel "Time [s]"
set ylabel "Postion [m]"
set title "Konato ODE Solver - Damped Harmonic Oscillator fv=2 k=100"
set xrange [ 0 : 2 ]
set yrange [ -1 : 1 ]
set mxtics 5
set mytics 5
set xtics 0.5
set ytics 0.5
plot "exdampfv2k100.dat" using 1:2 notitle w l
EOF

This would show the position x of the mass relative to the time t:

exdampfv2k100

Here is an other example of a graph produced by gnuplot for fv and k = 1:

dampedfv1k1

Conclusion

This tutorial showed the basic part of using Clojure and the com.konato.ode library to simulate ordinary differential equations. Examples provide with the library show extra features and usage of the library. For example, using a state space representation or a tranfer function is easy with the state2ode and transfer2ode functions. New Runge-Kunta method can be added. Fixed and Variable step are possibles. A analog control of a electric bicycle, a digital pid used with a system having non-linear caracteristics are shown and more. Subsequent tutorials will show more possibility of the library.

Usefull Clojure reference and thanks

Clojure main site

R. Mark Volkmann Insight

Stuart Halloway book Programming Clojure

Google discussion group on Clojure

The nice unit test framework test-is created by Stuart Sierra

, ,

5 Comments

com.konato.ode v1.00 is released

This is the first release of com.konato.ode an ordinary differentials equations solvers and simulation library in Clojure.
It does implements several fixed steps solvers: Euler, an Euler modified version, Runge-Kutta order 4 and a Runge-Kutta using parameters to implement RK4, Runge-Kutta-Fehlberg. It does implement a variable step ODE solver based on the popular rfk45 method. Others RK fixed or variable methods can be easily implemented using the RKP method. Functions are provided to convert State Space and Transfer functions to ODE who can be used by the system.

Current version is 1.00

Available at: git://github.com/konato/com.konato.ode.git

Tutorials at: Konato Blog

Methods based on: Gerald-Wheatley, 2003, Applied Numerical Analysis, Seventh Edition,
Pearson Education, ISBN 0-321-13304-8.

Copyright (c) 2009 Stephane Rousseau (stephaner gmail). All rights reserved.

Some simulations:

,

4 Comments