CDE

Task/Euler-method/0DESCRIPTION

@@ -0,0 +1,50 @@
+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 [[wp:Euler method|the wikipedia page]].
+The ODE has to be provided in the following form:
+
+:<math>\frac{dy(t)}{dt} = f(t,y(t))</math>
+
+with an initial value
+
+:<math>y(t_0) = y_0</math>
+
+To get a numeric solution, we replace the derivative on the LHS with a finite difference approximation:
+
+:<math>\frac{dy(t)}{dt}  \approx \frac{y(t+h)-y(t)}{h}</math>
+
+then solve for <math>y(t+h)</math>:
+
+:<math>y(t+h) \approx y(t) + h \, \frac{dy(t)}{dt}</math>
+
+which is the same as
+
+:<math>y(t+h) \approx y(t) + h \, f(t,y(t))</math>
+
+The iterative solution rule is then:
+
+:<math>y_{n+1} = y_n + h \, f(t_n, y_n)</math>
+
+<math>h</math> 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 <math>T(t_0) = T_0</math> cools down in an environment of temperature <math>T_R </math>:
+
+:<math>\frac{dT(t)}{dt} = -k \, \Delta T</math>
+
+or
+
+:<math>\frac{dT(t)}{dt} = -k \, (T(t) - T_R)</math>
+
+It says that the cooling rate <math>\frac{dT(t)}{dt}</math> of the object is proportional to the current temperature difference <math>\Delta T = (T(t) - T_R)</math> to the surrounding environment.
+
+The analytical solution, which we will compare to the numerical approximation, is
+
+:<math>T(t) = T_R + (T_0 - T_R) \; e^{-k t}</math>
+
+'''Task'''
+
+The task is to 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 2 s, 5 s and 10 s and to compare with the analytical solution.
+The initial temperature <math>T_0</math> shall be 100 °C, the room temperature <math>T_R</math> 20 °C, and the cooling constant <math>k</math> 0.07. The time interval to calculate shall be from 0 s to 100 s.
+
+A reference solution ([[#Common Lisp|Common Lisp]]) can be seen below. We see that bigger step sizes lead to reduced approximation accuracy.
+[[Image:Euler_Method_Newton_Cooling.png|center|750px]]

Task/Euler-method/1META.yaml

@@ -0,0 +1,2 @@
+---
+note: Mathematical operations

Task/Euler-method/ALGOL-68/euler-method.alg

@@ -0,0 +1,24 @@
+#
+Approximates y(t) in y'(t)=f(t,y) with y(a)=y0 and
+t=a..b and the step size h.
+#
+PROC euler = (PROC(REAL,REAL)REAL f, REAL y0, a, b, h)REAL: (
+    REAL y := y0,
+         t := a;
+    WHILE t < b DO
+      printf(($g(-6,3)": "g(-7,3)l$, t, y));
+      y +:= h * f(t, y);
+      t +:= h
+    OD;
+    printf($"done"l$);
+    y
+);
+
+# Example: Newton's cooling law #
+PROC newton cooling law = (REAL time, t)REAL: (
+    -0.07 * (t - 20)
+);
+
+main: (
+   euler(newton cooling law, 100, 0, 100,  10)
+)

Task/Euler-method/Ada/euler-method-1.ada

@@ -0,0 +1,11 @@
+generic
+   type Number is digits <>;
+package Euler is
+   type Waveform is array (Integer range <>) of Number;
+   function Solve
+            (  F      : not null access function (T, Y : Number) return Number;
+               Y0     : Number;
+               T0, T1 : Number;
+               N      : Positive
+            )  return Waveform;
+end Euler;

Task/Euler-method/Ada/euler-method-2.ada

@@ -0,0 +1,17 @@
+package body Euler is
+   function Solve
+            (  F      : not null access function (T, Y : Number) return Number;
+               Y0     : Number;
+               T0, T1 : Number;
+               N      : Positive
+            )  return Waveform is
+      dT : constant Number := (T1 - T0) / Number (N);
+   begin
+      return Y : Waveform (0..N) do
+         Y (0) := Y0;
+         for I in 1..Y'Last loop
+            Y (I) := Y (I - 1) + dT * F (T0 + dT * Number (I - 1), Y (I - 1));
+         end loop;
+      end return;
+   end Solve;
+end Euler;

Task/Euler-method/Ada/euler-method-3.ada

@@ -0,0 +1,18 @@
+with Ada.Text_IO;  use Ada.Text_IO;
+with Euler;
+
+procedure Test_Euler_Method is
+   package Float_Euler is new Euler (Float);
+   use Float_Euler;
+
+   function Newton_Cooling_Law (T, Y : Float) return Float is
+   begin
+      return -0.07 * (Y - 20.0);
+   end Newton_Cooling_Law;
+
+   Y : Waveform := Solve (Newton_Cooling_Law'Access, 100.0, 0.0, 100.0, 10);
+begin
+   for I in Y'Range loop
+      Put_Line (Integer'Image (10 * I) & ":" & Float'Image (Y (I)));
+   end loop;
+end Test_Euler_Method;

Task/Euler-method/BBC-BASIC/euler-method.bbc

@@ -0,0 +1,15 @@
+      PROCeuler("-0.07*(y-20)", 100, 0, 100, 2)
+      PROCeuler("-0.07*(y-20)", 100, 0, 100, 5)
+      PROCeuler("-0.07*(y-20)", 100, 0, 100, 10)
+      END
+
+      DEF PROCeuler(df$, y, a, b, s)
+      LOCAL t, @%
+      @% = &2030A
+      t = a
+      WHILE t <= b
+        PRINT t, y
+        y += s * EVAL(df$)
+        t += s
+      ENDWHILE
+      ENDPROC

Task/Euler-method/C++/euler-method.cpp

@@ -0,0 +1,32 @@
+#include <iomanip>
+#include <iostream>
+
+typedef double F(double,double);
+
+/*
+Approximates y(t) in y'(t)=f(t,y) with y(a)=y0 and
+t=a..b and the step size h.
+*/
+void euler(F f, double y0, double a, double b, double h)
+{
+    double y = y0;
+    for (double t = a; t < b; t += h)
+    {
+        std::cout << std::fixed << std::setprecision(3) << t << " " << y << "\n";
+        y += h * f(t, y);
+    }
+    std::cout << "done\n";
+}
+
+// Example: Newton's cooling law
+double newtonCoolingLaw(double, double t)
+{
+    return -0.07 * (t - 20);
+}
+
+int main()
+{
+    euler(newtonCoolingLaw, 100, 0, 100,  2);
+    euler(newtonCoolingLaw, 100, 0, 100,  5);
+    euler(newtonCoolingLaw, 100, 0, 100, 10);
+}

Task/Euler-method/C-sharp/euler-method.cs

@@ -0,0 +1,38 @@
+using System;
+
+namespace prog
+{
+	class MainClass
+	{
+		const float T0 = 100f;
+		const float TR = 20f;
+		const float k = 0.07f;
+		readonly static float[] delta_t = {2.0f,5.0f,10.0f};
+		const int n = 100;
+		
+		public delegate float func(float t);
+		static float NewtonCooling(float t)
+		{
+			return -k * (t-TR);			
+		}
+		
+		public static void Main (string[] args)
+		{
+			func f = new func(NewtonCooling);
+			for(int i=0; i<delta_t.Length; i++)
+			{
+				Console.WriteLine("delta_t = " + delta_t[i]);
+				Euler(f,T0,n,delta_t[i]);
+			}
+		}
+				
+		public static void Euler(func f, float y, int n, float h)
+		{
+			for(float x=0; x<=n; x+=h)
+			{
+				Console.WriteLine("\t" + x + "\t" + y);
+				y += h * f(y);	
+			}
+		}
+	}
+}

Task/Euler-method/C/euler-method-1.c

@@ -0,0 +1,44 @@
+#include <stdio.h>
+#include <math.h>
+
+typedef double (*deriv_f)(double, double);
+#define FMT " %7.3f"
+
+void ivp_euler(deriv_f f, double y, int step, int end_t)
+{
+	int t = 0;
+
+	printf(" Step %2d: ", (int)step);
+	do {
+		if (t % 10 == 0) printf(FMT, y);
+		y += step * f(t, y);
+	} while ((t += step) <= end_t);
+	printf("\n");
+}
+
+void analytic()
+{
+	double t;
+	printf("    Time: ");
+	for (t = 0; t <= 100; t += 10) printf(" %7g", t);
+	printf("\nAnalytic: ");
+
+	for (t = 0; t <= 100; t += 10)
+		printf(FMT, 20 + 80 * exp(-0.07 * t));
+	printf("\n");
+}
+
+double cooling(double t, double temp)
+{
+	return -0.07 * (temp - 20);
+}
+
+int main()
+{
+	analytic();
+	ivp_euler(cooling, 100, 2, 100);
+	ivp_euler(cooling, 100, 5, 100);
+	ivp_euler(cooling, 100, 10, 100);
+
+	return 0;
+}

Task/Euler-method/C/euler-method-2.c

@@ -0,0 +1,5 @@
+    Time:        0      10      20      30      40      50      60      70      80      90     100
+Analytic:  100.000  59.727  39.728  29.797  24.865  22.416  21.200  20.596  20.296  20.147  20.073
+ Step  2:  100.000  57.634  37.704  28.328  23.918  21.843  20.867  20.408  20.192  20.090  20.042
+ Step  5:  100.000  53.800  34.280  26.034  22.549  21.077  20.455  20.192  20.081  20.034  20.014
+ Step 10:  100.000  44.000  27.200  22.160  20.648  20.194  20.058  20.017  20.005  20.002  20.000

Task/Euler-method/Clay/euler-method.clay

@@ -0,0 +1,16 @@
+import printer.formatter as pf;
+
+euler(f, y, a, b, h) {
+    while (a < b) {
+        println(pf.rightAligned(2, a), " ", y);
+        a += h;
+        y += h * f(y);
+    }
+}
+
+main() {
+    for (i in [2.0, 5.0, 10.0]) {
+        println("\nFor delta = ", i, ":");
+        euler((temp) => -0.07 * (temp - 20), 100.0, 0.0, 100.0, i);
+    }
+}

Task/Euler-method/Common-Lisp/euler-method.lisp

@@ -0,0 +1,26 @@
+;; 't' usually means "true" in CL, but we need 't' here for time/temperature.
+(defconstant true 'cl:t)
+(shadow 't)
+
+
+;; Approximates y(t) in y'(t)=f(t,y) with y(a)=y0 and t=a..b and the step size h.
+(defun euler (f y0 a b h)
+
+  ;; Set the initial values and increments of the iteration variables.
+  (do ((t a  (incf t h))
+       (y y0 (incf y (* h (funcall f t y)))))
+
+      ;; End the iteration when t reaches the end b of the time interval.
+      ((>= t b) 'DONE)
+
+      ;; Print t and y(t) at every step of the do loop.
+      (format true "~6,3F  ~6,3F~%" t y)))
+
+
+;; Example: Newton's cooling law, f(t,T) = -0.07*(T-20)
+(defun newton-cooling (time T) (* -0.07 (- T 20)))
+
+;; Generate the data for all three step sizes (2,5 and 10).
+(euler #'newton-cooling 100 0 100  2)
+(euler #'newton-cooling 100 0 100  5)
+(euler #'newton-cooling 100 0 100 10)

Task/Euler-method/D/euler-method.d

@@ -0,0 +1,26 @@
+import std.stdio, std.range;
+
+/**
+Approximates y(t) in y'(t)=f(t,y) with y(a)=y0 and
+t=a..b and the step size h.
+*/
+void euler(F)(in F f, in double y0,
+              in double a, in double b, in double h) {
+    double y = y0;
+    foreach (t; iota(a, b, h)) {
+        writefln("%.3f  %.3f", t, y);
+        y += h * f(t, y);
+    }
+    writeln("done");
+}
+
+void main() {
+    /// Example: Newton's cooling law
+    static newtonCoolingLaw(in double time, in double t) {
+        return -0.07 * (t - 20);
+    }
+
+    euler(&newtonCoolingLaw, 100, 0, 100,  2);
+    euler(&newtonCoolingLaw, 100, 0, 100,  5);
+    euler(&newtonCoolingLaw, 100, 0, 100, 10);
+}

Task/Euler-method/Euler-Math-Toolbox/euler-method.euler

@@ -0,0 +1,21 @@
+>function dgleuler (f,x,y0) ...
+$  y=zeros(size(x)); y[1]=y0;
+$  for i=2 to cols(y);
+$  y[i]=y[i-1]+f(x[i-1],y[i-1])*(x[i]-x[i-1]);
+$  end;
+$  return y;
+$endfunction
+>function f(x,y) := -k*(y-TR)
+>k=0.07; TR=20; TS=100;
+>x=0:1:100; dgleuler("f",x,TS)[-1]
+ 20.0564137335
+>x=0:2:100; dgleuler("f",x,TS)[-1]
+ 20.0424631834
+>TR+(TS-TR)*exp(-k*TS)
+ 20.0729505572
+>x=0:5:100; plot2d(x,dgleuler("f",x,TS)); ...
+>  plot2d(x,TR+(TS-TR)*exp(-k*x),>add,color=red);
+>ode("f",x,TS)[-1] // Euler default solver LSODA
+ 20.0729505568
+>adaptiverunge("f",x,TS)[-1] // Adaptive Runge Method
+ 20.0729505572

Task/Euler-method/Fortran/euler-method.f

@@ -0,0 +1,57 @@
+program euler_method
+use iso_fortran_env, only: real64
+implicit none
+
+abstract interface
+  ! a derivative dy/dt as function of y and t
+  function derivative(y, t)
+    use iso_fortran_env, only: real64
+    real(real64) :: derivative
+    real(real64), intent(in) :: t, y
+  end function
+end interface
+
+real(real64), parameter :: T_0 = 100, T_room = 20, k = 0.07, a = 0, b = 100, &
+    h(3) = [2.0, 5.0, 10.0]
+
+integer :: i
+
+! loop over all step sizes
+do i = 1, 3
+  call euler(newton_cooling, T_0, a, b, h(i))
+end do
+
+contains
+
+! Approximates y(t) in y'(t) = f(y, t) with y(a) = y0 and t = a..b and the
+! step size h.
+subroutine euler(f, y0, a, b, h)
+  procedure(derivative) :: f
+  real(real64), intent(in) :: y0, a, b, h
+  real(real64) :: t, y
+
+  if (a > b) return
+  if (h <= 0) stop "negative step size"
+
+  print '("# h = ", F0.3)', h
+
+  y = y0
+  t = a
+
+  do
+    print *, t, y
+    t = t + h
+    if (t > b) return
+    y = y + h * f(y, t)
+  end do
+end subroutine
+
+
+! Example: Newton's cooling law, f(T, _) = -k*(T - T_room)
+function newton_cooling(T, unused) result(dTdt)
+  real(real64) :: dTdt
+  real(real64), intent(in) :: T, unused
+  dTdt = -k * (T - T_room)
+end function
+
+end program

Task/Euler-method/Go/euler-method.go

@@ -0,0 +1,68 @@
+package main
+
+import (
+    "fmt"
+    "math"
+)
+
+// fdy is a type for function f used in Euler's method.
+type fdy func(float64, float64) float64
+
+// eulerStep computes a single new value using Euler's method.
+// Note that step size h is a parameter, so a variable step size
+// could be used.
+func eulerStep(f fdy, x, y, h float64) float64 {
+    return y + h*f(x, y)
+}
+
+// Definition of cooling rate.  Note that this has general utility and
+// is not specific to use in Euler's method.
+
+// newCoolingRate returns a function that computes cooling rate
+// for a given cooling rate constant k.
+func newCoolingRate(k float64) func(float64) float64 {
+    return func(deltaTemp float64) float64 {
+        return -k * deltaTemp
+    }
+}
+
+// newTempFunc returns a function that computes the analytical solution
+// of cooling rate integrated over time.
+func newTempFunc(k, ambientTemp, initialTemp float64) func(float64) float64 {
+    return func(time float64) float64 {
+        return ambientTemp + (initialTemp-ambientTemp)*math.Exp(-k*time)
+    }
+}
+
+// newCoolingRateDy returns a function of the kind needed for Euler's method.
+// That is, a function representing dy(x, y(x)).
+//
+// Parameters to newCoolingRateDy are cooling constant k and ambient
+// temperature.
+func newCoolingRateDy(k, ambientTemp float64) fdy {
+    crf := newCoolingRate(k)
+    // note that result is dependent only on the object temperature.
+    // there are no additional dependencies on time, so the x parameter
+    // provided by eulerStep is unused.
+    return func(_, objectTemp float64) float64 {
+        return crf(objectTemp - ambientTemp)
+    }
+}
+
+func main() {
+    k := .07
+    tempRoom := 20.
+    tempObject := 100.
+    fcr := newCoolingRateDy(k, tempRoom)
+    analytic := newTempFunc(k, tempRoom, tempObject)
+    for _, deltaTime := range []float64{2, 5, 10} {
+        fmt.Printf("Step size = %.1f\n", deltaTime)
+        fmt.Println(" Time Euler's Analytic")
+        temp := tempObject
+        for time := 0.; time <= 100; time += deltaTime {
+            fmt.Printf("%5.1f %7.3f %7.3f\n", time, temp, analytic(time))
+            temp = eulerStep(fcr, time, temp, deltaTime)
+        }
+        fmt.Println()
+    }
+}

Task/Euler-method/Haskell/euler-method.hs

@@ -0,0 +1,15 @@
+import Text.Printf
+
+euler :: (Num a, Ord a) => (a -> a -> a) -> a -> a -> a -> a -> [(a,a)]
+euler f y0 a b h =
+	(a, y0) :
+	if a < b
+		then euler f (y0 + (f a y0) * h) (a + h) b h
+		else []
+
+newtonCooling :: Double -> Double -> Double
+newtonCooling _ t = -0.07 * (t - 20)
+
+main = do
+  mapM_ (uncurry $ printf "%6.3f %6.3f\n") $ euler newtonCooling 100 0 100 10
+  putStrLn "DONE"

Task/Euler-method/Java/euler-method.java

@@ -0,0 +1,33 @@
+public class Euler {
+  private static void euler (Callable f, double y0, int a, int b, int h) {
+    int t = a;
+    double y = y0;
+    while (t < b) {
+      System.out.println ("" + t + " " + y);
+      t += h;
+      y += h * f.compute (t, y);
+    }
+    System.out.println ("DONE");
+  }
+
+  public static void main (String[] args) {
+    Callable cooling = new Cooling ();
+    int[] steps = {2, 5, 10};
+    for (int stepSize : steps) {
+      System.out.println ("Step size: " + stepSize);
+      euler (cooling, 100.0, 0, 100, stepSize);
+    }
+  }
+}
+
+// interface used so we can plug in alternative functions to Euler
+interface Callable {
+  public double compute (int time, double t);
+}
+
+// class to implement the newton cooling equation
+class Cooling implements Callable {
+  public double compute (int time, double t) {
+    return -0.07 * (t - 20);
+  }
+}

Task/Euler-method/Lua/euler-method.lua

@@ -0,0 +1,22 @@
+T0 = 100
+TR = 20
+k = 0.07
+delta_t = { 2, 5, 10 }
+n = 100
+
+NewtonCooling = function( t ) return -k * ( t - TR ) end
+
+
+function Euler( f, y0, n, h )
+    local y = y0
+    for x = 0, n, h do
+	print( "", x, y )
+ 	y = y + h * f( y )
+    end
+end
+
+
+for i = 1, #delta_t do
+    print( "delta_t = ", delta_t[i] )
+    Euler( NewtonCooling, T0, n, delta_t[i] )
+end

Task/Euler-method/Perl/euler-method.pl

@@ -0,0 +1,32 @@
+sub euler_method {
+        my ($t0, $t1, $k, $step_size) = @_;
+        my @results = ( [0, $t0] );
+
+        for (my $s = $step_size; $s <= 100; $s += $step_size) {
+                $t0 -= ($t0 - $t1) * $k * $step_size;
+                push @results, [$s, $t0];
+        }
+
+        return @results;
+}
+
+sub analytical {
+        my ($t0, $t1, $k, $time) = @_;
+        return ($t0 - $t1) * exp(-$time * $k) + $t1
+}
+
+my ($T0, $T1, $k) = (100, 20, .07);
+my @r2  = grep { $_->[0] % 10 == 0 } euler_method($T0, $T1, $k, 2);
+my @r5  = grep { $_->[0] % 10 == 0 } euler_method($T0, $T1, $k, 5);
+my @r10 = grep { $_->[0] % 10 == 0 } euler_method($T0, $T1, $k, 10);
+
+print "Time\t      2     err(%)      5     err(%)    10      err(%)  Analytic\n", "-" x 76, "\n";
+for (0 .. $#r2) {
+        my $an = analytical($T0, $T1, $k, $r2[$_][0]);
+        printf "%4d\t".("%9.3f" x 7)."\n",
+                $r2 [$_][0],
+                $r2 [$_][1], ($r2 [$_][1] / $an) * 100 - 100,
+                $r5 [$_][1], ($r5 [$_][1] / $an) * 100 - 100,
+                $r10[$_][1], ($r10[$_][1] / $an) * 100 - 100,
+                $an;
+}

Task/Euler-method/PicoLisp/euler-method.l

@@ -0,0 +1,14 @@
+(load "@lib/math.l")
+
+(de euler (F Y A B H)
+   (while (> B A)
+      (prinl (round A) " " (round Y))
+      (inc 'Y (*/ H (F A Y) 1.0))
+      (inc 'A H) ) )
+
+(de newtonCoolingLaw (A B)
+   (*/ -0.07 (- B 20.) 1.0) )
+
+(euler newtonCoolingLaw 100.0 0 100.0 2.0)
+(euler newtonCoolingLaw 100.0 0 100.0 5.0)
+(euler newtonCoolingLaw 100.0 0 100.0 10.0)

Task/Euler-method/Python/euler-method.py

@@ -0,0 +1,11 @@
+def euler(f,y0,a,b,h):
+	t,y = a,y0
+	while t < b:
+		print "%6.3f %6.3f" % (t,y)
+		t += h
+		y += h * f(t,y)
+
+def newtoncooling(time, temp):
+	return -0.07 * (temp - 20)
+
+euler(newtoncooling,100,0,100,10)

Task/Euler-method/Ruby/euler-method.rb

@@ -0,0 +1,10 @@
+def euler(y0,a,b,h, &block)
+  t,y = a,y0
+  while t < b
+    puts "%6.3f %6.3f" % [t,y]
+    t += h
+    y += h * block.call(t,y)
+  end
+end
+
+euler(100,0,100,10) {|time, temp| -0.07 * (temp - 20) }

Task/Euler-method/Scala/euler-method.scala

@@ -0,0 +1,29 @@
+object App{
+
+  def main(args : Array[String]) = {
+
+    def cooling( step : Int ) = {
+      eulerStep( (step , y) => {-0.07 * (y - 20)} ,
+        100.0,0,100,step)
+    }
+    cooling(10)
+    cooling(5)
+    cooling(2)
+  }
+  def eulerStep( func : (Int,Double) => Double,y0 : Double,
+    begin : Int, end : Int , step : Int) = {
+
+    println("Step size: %s".format(step))
+
+    var current : Int = begin
+    var y : Double = y0
+    while( current <= end){
+      println( "%d %.5f".format(current,y))
+      current += step
+      y += step * func(current,y)
+    }
+
+    println("DONE")
+  }
+
+}

Task/Euler-method/Tcl/euler-method-1.tcl

@@ -0,0 +1,9 @@
+proc euler {f y0 a b h} {
+    puts "computing $f over \[$a..$b\], step $h"
+    set y [expr {double($y0)}]
+    for {set t [expr {double($a)}]} {$t < $b} {set t [expr {$t + $h}]} {
+	puts [format "%.3f\t%.3f" $t $y]
+	set y [expr {$y + $h * double([$f $t $y])}]
+    }
+    puts "done"
+}

Task/Euler-method/Tcl/euler-method-2.tcl

@@ -0,0 +1,7 @@
+proc newtonCoolingLaw {time temp} {
+    expr {-0.07 * ($temp - 20)}
+}
+
+euler newtonCoolingLaw 100 0 100 2
+euler newtonCoolingLaw 100 0 100 5
+euler newtonCoolingLaw 100 0 100 10