Mercurial > repos > public > sbplib

comparison +time/Test.m @ 1114:f2988a63c3aa feature/timesteppers

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add tests for time.Rungekutta4 and time.Rungekutta4SecondOrder
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 10 Apr 2019 23:10:13 +0200
parents
children
comparison
equal deleted inserted replaced
1113:47e86b5270ad 1114:f2988a63c3aa
1 function tests = Test()
2 tests = functiontests(localfunctions);
3 end
4
5 function testRungekutta4(tc)
6 F = @(t,y) -0.5*y;
7 dt = 0.01;
8 t0 = 0;
9 y0 = 1;
10 T = 2;
11
12 ts1 = time.Rungekutta4(F, dt, t0, y0);
13 ts1.stepN(T/dt);
14
15 ts2 = time.Rungekutta4(F, dt/2, t0, y0);
16 ts2.stepN(T/(dt/2));
17
18 y_true = exp(-0.5*T);
19
20 e1 = y_true - ts1.getV();
21 e2 = y_true - ts2.getV();
22
23 tc.verifyEqual(e1/e2, 2^4, 'AbsTol', 0.5);
24 end
25
26 function testRungekutta4SecondOrder(tc)
27 F = @(t,y,y_t) -4*y - 1/4*y_t;
28 dt = 0.01;
29 t0 = 0;
30 y0 = 1;
31 y0t = 1/2;
32 T = 2;
33
34 ts1 = time.Rungekutta4SecondOrder(F, dt, t0, y0, y0t);
35 ts1.stepN(T/dt);
36
37 ts2 = time.Rungekutta4SecondOrder(F, dt/2, t0, y0, y0t);
38 ts2.stepN(T/(dt/2));
39
40 y_true = 1/51*exp(-T/8)*(sqrt(255)*sin(sqrt(255)/8*T) + 51*cos(sqrt(255)/8*T));
41
42 e1 = y_true - ts1.getV();
43 e2 = y_true - ts2.getV();
44
45 tc.verifyEqual(e1/e2, 2^4, 'AbsTol', 0.5);
46 end