Euler method: Difference between revisions

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]]

 

=={{header|Ada}}==

The solution is generic, usable for any floating point type. The package specification:

generic

type Number is digits <>;

) return Waveform;

end Euler;

The function Solve returns the solution of the differential equation for each of N+1 points, starting from the point T0. The implementation:

package body Euler is

function Solve

end Solve;

end Euler;

The test program:

with Ada.Text_IO; use Ada.Text_IO;

with Euler;

end loop;

end Test_Euler_Method;

Sample output:

<pre>

100: 2.00005E+01

</pre>

=={{header|ALGOL 68}}==

{{trans|D}} Note: This specimen retains the original [[#D|D]] coding style.

{{works with|ALGOL 68G|Any - tested with release [http://sourceforge.net/projects/algol68/files/algol68g/algol68g-1.18.0/algol68g-1.18.0-9h.tiny.el5.centos.fc11.i386.rpm/download 1.18.0-9h.tiny].}}

{{wont work with|ELLA ALGOL 68|Any (with appropriate job cards) - tested with release [http://sourceforge.net/projects/algol68/files/algol68toc/algol68toc-1.8.8d/algol68toc-1.8-8d.fc9.i386.rpm/download 1.8-8d] - due to extensive use of '''format'''[ted] ''transput''.}}

Approximates y(t) in y'(t)=f(t,y) with y(a)=y0 and

t=a..b and the step size h.

main: (

euler(newton cooling law, 100, 0, 100, 10)

Ouput:

<pre>

done

</pre>

 

=={{header|BASIC}}==

==={{header|BBC BASIC}}===

PROCeuler("-0.07*(y-20)", 100, 0, 100, 5)

PROCeuler("-0.07*(y-20)", 100, 0, 100, 10)

t += s

ENDWHILE

'''Output:'''

<pre>

100.000 20.000

</pre>

 

==={{header|FreeBASIC}}===

'Custom rounding

#define round(x,N) Rtrim(Rtrim(Left(Str((x)+(.5*Sgn((x)))/(10^(N))),Instr(Str((x)+(.5*Sgn((x)))/(10^(N))),".")+(N)),"0"),".")

Euler(-.07*(y-20),100,0,100,10,"print")

Sleep

outputs (steps 5 and 10)

<pre>

 

</pre>

 

==={{header|Run BASIC}}===

x = euler-0.07,-20, 100, 0, 100, 5)

x = euler(-0.07,-20, 100, 0, 100, 10)

t = t + s

WEND

<pre>===== da:-0.07 db:-20 y:100 a:0 b:100 s:2 ===================

0 100

70 20.017496

80 20.0052488</pre>

 

=={{header|C}}==

#include <math.h>

 

 

return 0;

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

 

=={{header|C++}}==

{{trans|D}}

#include <iostream>

 

euler(newtonCoolingLaw, 100, 0, 100, 5);

euler(newtonCoolingLaw, 100, 0, 100, 10);

Last part of output:

<pre>

done

</pre>

 

=={{header|Clay}}==

 

import printer.formatter as pf;

 

}

}

 

Example output:

80 20.0052488

90 20.00157464

</pre>

 

{{works with|Visual COBOL}}

The following is in the Managed COBOL dialect:

PROCEDURE DIVISION USING VALUE t AS FLOAT-LONG

RETURNING ret AS FLOAT-LONG.

END-PERFORM

END METHOD.

 

Example output:

x = 80.0000, y = 20.0052

x = 90.0000, y = 20.0016

</pre>

 

=={{header|Common Lisp}}==

(defconstant true 'cl:t)

(shadow 't)

(euler #'newton-cooling 100 0 100 2)

(euler #'newton-cooling 100 0 100 5)

 

 

(defun newton-cooling (time temperature)

(loop for time from a below b by h

for y = y0 then (+ y (* h (funcall f time y)))

 

<pre>

 

=={{header|D}}==

 

/// Approximates y(t) in y'(t)=f(t,y) with y(a)=y0 and t=a..b and the step size h.

euler(newtonCoolingLaw, 100, 0, 100, 5);

euler(newtonCoolingLaw, 100, 0, 100, 10);

Last part of the output:

<pre>...

90.000 20.002

done</pre>

 

=={{header|Delphi}}==

 

[[Euler_method#Pascal | Pascal]]

 

=={{header|Elixir}}==

{{trans|Ruby}}

def method(_, _, t, b, _) when t>b, do: :ok

def method(f, y, t, b, h) do

IO.puts "\nStep = #{step}"

Euler.method(f, 100.0, 0.0, 100.0, step)

 

{{out}}

=={{header|Erlang}}==

 

-module(euler).

-export([main/0, euler/5]).

euler(fun cooling/2, 100, 0, 10, 100),

ok.

{{out}}

<pre>

=={{header|Euler Math Toolbox}}==

 

>function dgleuler (f,x,y0) ...

$ y=zeros(size(x)); y[1]=y0;

>adaptiverunge("f",x,TS)[-1] // Adaptive Runge Method

20.0729505572

 

=={{header|F_Sharp|F#}}==

(t0, y0)

|> Seq.unfold (fun (t, y) -> Some((t,y), ((t + h), (y + h * (f t y)))))

|> Seq.takeWhile (fun (t,_) -> t <= b)

|> Seq.iter (printfn "%A")

Output for the above (step size 10)

<pre>(0.0, 100.0)

(90.0, 20.00157464)

(100.0, 20.00047239)</pre>

 

=={{header|Forth}}==

20e f- -0.07e f* ;

 

100e ' newton-cooling-law 2 100 euler cr

100e ' newton-cooling-law 5 100 euler cr

 

=={{header|Fortran}}==

{{works with|Fortran|2008}}

use iso_fortran_env, only: real64

implicit none

end function

 

Output for <code>h = 10</code>:

<pre>

determined temperature (so we are computing the error).

 

20.0 + (t0 - 20.0) * f64.exp(-0.07*time)

 

-0.07 * (temperature-20.0)

 

let temps = replicate steps 0.0

let temps[i] = f64.abs(t-analytic t0 x)

in (t + h * cooling x t,

temps)

in temps

 

=={{header|Go}}==

 

import (

fmt.Println()

}

Output, truncated:

<pre>

 

The ''eulerMapping'' closure produces an entire approximating sequence, expressed as a Map object. Here, ''x₀'' and ''y₀'' together are the first point in the sequence, the ODE initial conditions. ''h'' and ''dydx'' are again the step distance and the derivative closure. ''stopCond'' is a closure representing a "stop condition" that causes the the ''eulerMapping'' closure to stop after a finite number of steps; the given default value causes ''eulerMapping'' to stop after 100 steps.

(yn + h * dydx(xn, yn)) as BigDecimal

}

yMap

}

 

 

'''Specific Euler Method Solution for the "Temperature Diffusion" Problem''' (with Newton's derivative formula and constants for environment temperature and object conductivity given)

assert dtdsNewton.maximumNumberOfParameters == 4

 

 

def tEulerH = eulerMapping.rcurry(dtds) { s, t -> s >= 100 }

 

 

'''Newton's Analytic Temperature Diffusion Solution''' (for comparison)

tR + (t0 - tR) * Math.exp(k * (s0 - s))

}

 

def tAnalytic = tNewton.rcurry(0, 100, 20, 0.07)

 

 

'''Specific solutions for 3 step sizes''' (and initial time and temperature)

def tEuler = tEulerH.rcurry(h)

assert tEuler.maximumNumberOfParameters == 2

printf('%5.0f %8.4f %8.4f %9.6f\n', s, tA, tE, relError)

}

 

 

=={{header|Haskell}}==

Modular solution which separates the solver and a method. Moreover it works on a given mesh which can be irregular.

dsolveBy _ _ [] _ = error "empty solution interval"

dsolveBy method f mesh x0 = zip mesh results

where results = scanl (method f) x0 intervals

 

It is better to use strict <code>Data.List.scanl'</code> in the solver but avoiding highlighting problems we leave lazy <code>scanl</code> function.

Some possible methods:

 

euler f x (t1,t2) = x + (t2 - t1) * f t1 x

 

k3 = f (t1 + h/2) (x + h/2*k2)

k4 = f (t1 + h) (x + h*k3)

 

Graphical output, using EasyPlot:

 

 

newton t temp = -0.07 * (temp - 20)

where sol1 = dsolveBy euler newton [0,10..100] 100

sol2 = dsolveBy rk2 newton [0,10..100] 100

 

=={{header|Icon}} and {{header|Unicon}}==

This solution works in both Icon and Unicon. It takes advantage of the <code>proc</code> procedure, which converts a string naming a procedure into a call to that procedure.

 

invocable "newton_cooling" # needed to use the 'proc' procedure

 

euler ("newton_cooling", 100, 0, 100, step_size)

end

 

Sample output:

=={{header|J}}==

'''Solution:'''

euler=: adverb define

'Y0 a b h'=. 4{. y

 

ncl=: _0.07 * -&20 NB. Newton's Cooling Law

'''Example:'''

... NB. output redacted for brevity

ncl euler 100 0 100 5

80 20.0052

90 20.0016

 

=={{header|Java}}==

 

public class Euler {

private static void euler (Callable f, double y0, int a, int b, int h) {

}

}

 

Output for step = 10;

=={{header|JavaScript}}==

{{trans|Python}}

// Function that takes differential-equation, initial condition,

// ending x, and step size as parameters

 

eulersMethod(cooling, 0, 100, 100, 10);

=={{header|jq}}==

{{works with|jq|1.4}}

# (x1,y1), ending x (x2), and step size as parameters;

# it emits the y values at each iteration.

# The following solves the task:

# (2,5,10) | [., [ euler_method(cooling; 0; 100; 100; .) ] ]

For brevity, we modify euler_method so that it only shows the final value of y:

def euler_solution(df; x1; y1; x2; h):

def recursion(exp): reduce recurse(exp) as $x (.; $x);

else empty

end )

'''Example''':

{{Out}}

[1,[20.05641373347389]]

[2,[20.0424631833732]]

[5,[20.01449963666907]]

[10,[20.000472392]]

 

=={{header|Julia}}==

 

# Prints a series of arbitrary values in a tabular form, left aligned in cells with a given width

t = 0

for T in euler(y -> -0.07 * (y - 20.0), 100.0, 0, 100, h)

t += h

end

println()

 

{{out}}

=={{header|Kotlin}}==

{{trans|C}}

 

typealias Deriv = (Double) -> Double // only one parameter needed here

for (i in listOf(2, 5, 10))

euler(::cooling, 100.0, i, 100)

 

{{out}}

Step 10: 100.000 44.000 27.200 22.160 20.648 20.194 20.058 20.017 20.005 20.002 20.000

</pre>

 

=={{header|Lua}}==

TR = 20

k = 0.07

Euler( NewtonCooling, T0, n, delta_t[i] )

end

 

=={{header|Mathematica}} / {{header|Wolfram Language}}==

Better methods for differential equation solving are built into Mathematica, so the typical user would omit the Method and StartingStepSize options in the code below. However since the task requests Eulers method, here is the bad solution...

{{out}}

<pre>euler[2, 100]

euler[10, 100]

20.0005</pre>

 

=={{header|Maxima}}==

[t: a, y: y0, tg: [a], yg: [y0]],

unless t>=b do (

[xlabel, "t / seconds"],

[ylabel, "Temperature / C"]);

 

=={{header|МК-61/52}}==

+ П4 С/П ИП2 ИП3 + П2 БП 05 ...

 

Instead of dots typed calculation program equation ''f(u, t)'', where the arguments are ''t'' = Р2, ''u'' = Р4.

 

=={{header|Nim}}==

 

proc euler(f: proc (x,y: float): float; y0, a, b, h: float) =

y += h * f(t,y)

 

-0.07 * (temp - 20)

 

<pre>0.000 100.000

10.000 44.000

 

=={{header|Objeck}}==

class EulerMethod {

T0 : static : Float;

}

}

 

Output:

40 20.648

</pre>

 

=={{header|OCaml}}==

 

* Given a function, and stepsize, provides a function of (t,y) which

* returns the next step: (t',y'). *)

 

(* analytic solution for Newton cooling *)

 

Using the above functions to produce the task results:

let cool = euler (newton_cooling ~k:0.07 ~tr:20.)

 

results (cool ~step:10.)

results (cool ~step:5.)

 

Example output:

=={{header|Oforth}}==

 

| t |

a b h step: t [

System.Out t <<wjp(6, JUSTIFY_RIGHT, 3) " : " << y << cr

t y f perform h * y + ->y

 

Usage :

 

y 20 - -0.07 * ;

 

euler(#newtonCoolingLaw, 100.0, 0.0, 100.0, 2)

euler(#newtonCoolingLaw, 100.0, 0.0, 100.0, 5)

 

{{out}}

100 : 20.000472392

</pre>

 

=={{header|Pascal}}==

Euler code for Free Pascal - Delphi mode. Apart from the function-pointer calling convention for the NewtonCooling method, this example is ISO-7185 standard Pascal.

 

 

{$mode delphi}

END.

 

 

 

 

=={{header|Perl}}==

my ($t0, $t1, $k, $step_size) = @_;

my @results = ( [0, $t0] );

$an;

}

Output:<pre>Time 2 err(%) 5 err(%) 10 err(%) Analytic

----------------------------------------------------------------------------

90 20.090 -0.281 20.034 -0.559 20.002 -0.721 20.147

100 20.042 -0.152 20.014 -0.291 20.000 -0.361 20.073</pre>

 

=={{header|Phix}}==

{{out}}

<pre>

 

=={{header|PicoLisp}}==

 

(de euler (F Y A B H)

(euler newtonCoolingLaw 100.0 0 100.0 2.0)

(euler newtonCoolingLaw 100.0 0 100.0 5.0)

Output:

<pre>...

 

=={{header|PL/I}}==

 

declare (x, y, z) float;

end f;

 

 

Only the final two outputs are shown, for brevity.

100 20.0729503632 20.0004730225

</pre>

 

=={{header|PowerShell}}==

{{works with|PowerShell|4.0}}

function euler (${f}, ${y}, $y0, $t0, $tEnd) {

function f-euler ($tn, $yn, $h) {

}

euler f y $y0 $t0 $tEnd | Format-Table -AutoSize

<b>Output:</b>

<pre>

 

=={{header|PureBasic}}==

Prototype.d Func(Time, t)

 

Print(#CRLF$ + #CRLF$ + "Press ENTER to exit"): Input()

CloseConsole()

<pre>...

85.000 20.053

=={{header|Python}}==

{{trans|Common Lisp}}

t,y = a,y0

while t <= b:

 

euler(newtoncooling,100,0,100,10)

Output:

<pre>

80.000 20.005

90.000 20.002

</pre>

 

=={{header|R}}==

{{trans|Python}}

{

t <- a

}

 

return(-0.07*(temp-20))

 

Output:

<pre> 0.000 100.000

 

The ODE solver:

(define (ODE-solve f init

#:x-max x-max

(reverse

(iterate-while (λ (x . y) (<= x x-max)) (method f h) init)))

 

It uses the default integration method <tt>euler</tt>, defined separately.

 

(define (euler F h)

(λ (x y) (list (+ x h) (+ y (* h (F x y))))))

 

A general-purpose procedure which evalutes a given function ''f'' repeatedly starting with argument ''x'', while all results satisfy a predicate ''test''. Returns a list of iterations.

 

(define (iterate-while test f x)

(let next ([result x]

(next (apply f result) (cons result list-of-results))

list-of-results)))

 

Textual output:

> (define (newton-cooling t T)

(* -0.07 (- T 20)))

(90 20.00157464)

(100 20.000472392))

 

Plotting results:

> (require plot)

> (plot

'(red blue black))

#:legend-anchor 'top-right)

[[File:euler1.jpg]]

 

[http://en.wikipedia.org/wiki/Midpoint_method 2-nd order Runge-Kutta method]:

 

(define (RK2 F h)

(λ (x y)

(list (+ x h) (+ y (* h (F (+ x (* 1/2 h))

(+ y (* 1/2 h (F x y)))))))))

 

[http://en.wikipedia.org/wiki/Adams_method#Two-step_Adams.E2.80.93Bashforth Two-step Adams–Bashforth method]

(define (adams F h)

(case-lambda

(let ([f (F x y)])

(list (+ x h) (+ y (* 3/2 h f) (* -1/2 h f′)) f))]))

 

[http://en.wikipedia.org/wiki/Adaptive_stepsize Adaptive one-step method] modifier using absolute accuracy ''ε''

(define ((adaptive method ε) F h0)

(case-lambda

(list x1 (+ y1 τ) (* 2 h′)))]))

 

 

Comparison of different integration methods

> (define (solve-newton-cooling-by m)

(ODE-solve newton-cooling '(0 100)

#:color 'blue #:label "Adaptive Runge-Kutta"))

#:legend-anchor 'top-right)

 

[[File:euler2.jpg]]

 

See also [[Runge-Kutta method#Racket]]

 

=={{header|REXX}}==

===version 1===

{{trans|PLI}}

* 24.05.2013 Walter Pachl translated from PL/I

**********************************************************************/

r=r+0

Return r

Output:

<pre>

===version 2===

This REXX version allows values to be specified via the command line (CL).

<br>It also shows the percentage difference (analytic vs. Euler's method) for each calculation.

parse arg Ti Tr cc tt ss /*obtain optional arguments from the CL*/

@= '═' /*the character used in title separator*/

exit /*stick a fork in it, we're all done. */

/*──────────────────────────────────────────────────────────────────────────────────────*/

{{out|output|text=&nbsp; when using the default inputs:}}

<pre>

time in

══seconds══ ══Euler method══ ═════analytic═════ ══difference══

 

 

time in

══seconds══ ══Euler method══ ═════analytic═════ ══difference══

 

 

time in

══seconds══ ══Euler method══ ═════analytic═════ ══difference══

</pre>

 

=={{header|Ring}}==

decimals(3)

see euler("return -0.07*(y-20)", 100, 0, 100, 2) + nl

t += s

end

Output:

<pre>

10 57.634

</pre>

 

=={{header|Ruby}}==

{{trans|Python}}

a.step(b,h) do |t|

puts "%7.3f %7.3f" % [t,y]

euler(100,0,100,step) {|time, temp| -0.07 * (temp - 20) }

puts

{{out}}

<pre style="height: 40ex; overflow: scroll">

=={{header|Rust}}==

{{trans|Kotlin}}

print!(" Time: ");

for t in (0..100).step_by(10) {

euler(|temp| -0.07 * (temp - 20.0), 100.0, i, 100);

}

{{out}}

<pre> Time: 0 10 20 30 40 50 60 70 80 90

=={{header|Scala}}==

 

object App{

 

 

}

 

Output for step = 10;

</pre>

 

 

=={{header|SequenceL}}==

 

import <Utilities/Sequence.sl>;

 

output when x > n

else

Based on C# version [http://rosettacode.org/wiki/Euler_method#C.23] but using tail recursion instead of looping.

{{out}}

=={{header|Sidef}}==

{{trans|Perl}}

var results = [[0, t0]]

for s in (step_size..100 -> by(step_size)) {

r10[i][1], (r10[i][1] / an) * 100 - 100,

an)

{{out}}

<pre>

 

=={{header|Smalltalk}}==

| t y |

t := a.

 

ODESolver new eulerOf: [:time :temp| -0.07 * (temp - 20)] init: 100 from: 0 to: 100 step: 10

Transcript:

<pre>

80 20.005

90 20.002

</pre>

 

=={{header|Tcl}}==

{{trans|C++}}

puts "computing $f over \[$a..$b\], step $h"

set y [expr {double($y0)}]

}

puts "done"

Demonstration with the Newton Cooling Law:

expr {-0.07 * ($temp - 20)}

}

euler newtonCoolingLaw 100 0 100 2

euler newtonCoolingLaw 100 0 100 5

End of output:

<pre>

90.000 20.002

done

</pre>

 

=={{header|VBA}}==

Dim t As Integer

Debug.Print " Step "; step; ": ",

ivp_euler r_cooling, 100, 5, 100

ivp_euler r_cooling, 100, 10, 100

<pre> 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 5 : 100,000 53,800 34,281 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 </pre>

=={{header|XPL0}}==

 

proc Euler(Step); \Display cooling temperatures using Euler's method

 

Format(3,2); \display cooling temps using differential eqn.

Time:= 0.0;

while Time <= 100.1 do

Euler(5);

Euler(10);

 

Output:

Step 5 100.00 53.80 34.28 26.03 22.55 21.08 20.46 20.19 20.08 20.03 20.01

Step 10 100.00 44.00 27.20 22.16 20.65 20.19 20.06 20.02 20.01 20.00 20.00

</pre>

 

=={{header|zkl}}==

{{trans|C}}

fcn ivp_euler(f,y,step,end_t){

ivp_euler(cooling, 100.0, 2, 100);

ivp_euler(cooling, 100.0, 5, 100);

{{out}}

<pre>

=={{header|ZX Spectrum Basic}}==

{{trans|BBC_BASIC}}

20 LET t=a