Problem Statement

Euler's method numerically approximates solutions of first-order ordinary differential equations (ODEs) with a given initial value.   It is an explicit method for solving initial value problems (IVPs), as described in the wikipedia page. The ODE has to be provided in the following form: with an initial value To get a numeric solution, we replace the derivative on the   LHS   with a finite difference approximation: then solve for

y ( t + h )

{\displaystyle y(t+h)}

: which is the same as The iterative solution rule is then: where

h

{\displaystyle h}

is the step size, the most relevant parameter for accuracy of the solution.   A smaller step size increases accuracy but also the computation cost, so it has always has to be hand-picked according to the problem at hand.

Example: Newton's Cooling Law Newton's cooling law describes how an object of initial temperature

T (

t

0

)

T

0

{\displaystyle T(t_{0})=T_{0}}

cools down in an environment of temperature

T

R

{\displaystyle T_{R}}

: or

It says that the cooling rate

d T ( t )

d t

{\displaystyle {\frac {dT(t)}{dt}}}

of the object is proportional to the current temperature difference

Δ T

( T ( t ) −

T

R

)

{\displaystyle \Delta T=(T(t)-T_{R})}

to the surrounding environment. The analytical solution, which we will compare to the numerical approximation, is

Implement a routine of Euler's method and then to use it to solve the given example of Newton's cooling law with it for three different step sizes of: and to compare with the analytical solution.

A reference solution (Common Lisp) can be seen below.   We see that bigger step sizes lead to reduced approximation accuracy.

Let's start with the solution:

(ns newton-cooling
  (:gen-class))

(defn euler [f y0 a b h]
  "Euler's Method.
  Approximates y(time) in y'(time)=f(time,y) with y(a)=y0 and t=a..b and the step size h."
  (loop [t a
         y y0
         result []]
    (if (<= t b)
        (recur (+ t h) (+ y (* (f (+ t h) y) h)) (conj result [(double t) (double y)]))
        result)))

(defn newton-coolling [t temp]
  "Newton's cooling law, f(t,T) = -0.07*(T-20)"
  (* -0.07 (- temp 20)))

; Run for case h = 10
(println "Example output")
(doseq [q (euler newton-coolling 100 0 100 10)]
  (println (apply format "%.3f %.3f" q)))