Euler method

Euler method
源自於 https://rosettacode.org/wiki/Euler_method#C.2B.2B

You are encouraged to solve this taskaccording to the task description, using any language you may know.

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:

${\displaystyle {\frac {dy(t)}{dt}}=f(t,y(t))}$

with an initial value

${\displaystyle y(t_{0})=y_{0}}$

To get a numeric solution, we replace the derivative on the   LHS   with a finite difference approximation:

${\displaystyle {\frac {dy(t)}{dt}}\approx {\frac {y(t+h)-y(t)}{h}}}$

then solve for ${\displaystyle y(t+h)}$:

${\displaystyle y(t+h)\approx y(t)+h\,{\frac {dy(t)}{dt}}}$

which is the same as

${\displaystyle y(t+h)\approx y(t)+h\,f(t,y(t))}$

The iterative solution rule is then:

${\displaystyle y_{n+1}=y_{n}+h\,f(t_{n},y_{n})}$

where   ${\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   ${\displaystyle T(t_{0})=T_{0}}$   cools down in an environment of temperature   ${\displaystyle T_{R}}$:

${\displaystyle {\frac {dT(t)}{dt}}=-k\,\Delta T}$

or

${\displaystyle {\frac {dT(t)}{dt}}=-k\,(T(t)-T_{R})}$

It says that the cooling rate   ${\displaystyle {\frac {dT(t)}{dt}}}$   of the object is proportional to the current temperature difference   ${\displaystyle \Delta T=(T(t)-T_{R})}$   to the surrounding environment.

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

${\displaystyle T(t)=T_{R}+(T_{0}-T_{R})\;e^{-kt}}$

Task

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.

Initial values

•   initial temperature   ${\displaystyle T_{0}}$   shall be   100 °C
•   room temperature   ${\displaystyle T_{R}}$   shall be   20 °C
•   cooling constant     ${\displaystyle k}$     shall be   0.07
•   time interval to calculate shall be from   0 s   ──►   100 s

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

C++程式

#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);

}

輸出結果

0.000 100.000

2.000 88.800

4.000 79.168

6.000 70.884

8.000 63.761

10.000 57.634

12.000 52.365

14.000 47.834

16.000 43.937

18.000 40.586

20.000 37.704

22.000 35.226

24.000 33.094

26.000 31.261

28.000 29.684

30.000 28.328

32.000 27.163

34.000 26.160

36.000 25.297

38.000 24.556

40.000 23.918

42.000 23.369

44.000 22.898

46.000 22.492

48.000 22.143

50.000 21.843

52.000 21.585

54.000 21.363

56.000 21.172

58.000 21.008

60.000 20.867

62.000 20.746

64.000 20.641

66.000 20.551

68.000 20.474

70.000 20.408

72.000 20.351

74.000 20.302

76.000 20.259

78.000 20.223

80.000 20.192

82.000 20.165

84.000 20.142

86.000 20.122

88.000 20.105

90.000 20.090

92.000 20.078

94.000 20.067

96.000 20.057

98.000 20.049

done

0.000 100.000

5.000 72.000

10.000 53.800

15.000 41.970

20.000 34.280

25.000 29.282

30.000 26.034

35.000 23.922

40.000 22.549

45.000 21.657

50.000 21.077

55.000 20.700

60.000 20.455

65.000 20.296

70.000 20.192

75.000 20.125

80.000 20.081

85.000 20.053

90.000 20.034

95.000 20.022

done

0.000 100.000

10.000 44.000

20.000 27.200

30.000 22.160

40.000 20.648

50.000 20.194

60.000 20.058

70.000 20.017

80.000 20.005

90.000 20.002

done