Skip to content

2.2. Ordinary differential equations

Let the source function \(F:[t_0, T] \times ℝ \to ℝ\) be given, and the initial data \(u_0:ℝ \to ℝ\) be given. Consider the initial value problem:

\[ u'(t) = F(t, u(t)) \]

with the initial condition \(u(t_0) = u_0(t_0)\).

To solve the problem numerically, the subpackage specular.ode provides the following numerical schemes:

  • the specular Euler scheme (Type 1 ~ 6)
  • the specular trigonometric scheme
  • the specular Heun scheme
  • the explicit Euler scheme
  • the implicit Euler scheme
  • the Crank-Nicolson scheme

2.2.1. Specular Euler scheme

All functions return an instance of the ODEResult class that encapsulates the numerical results.

import specular

def F(t, u):
    return -2*u 

specular.Euler_scheme(of_Type='1', F=F, t_0=0.0, u_0=1.0, T=2.5, h=0.1)

Running the specular Euler scheme of Type 1: 100%|██████████| 24/24 [00:00<?, ?it/s]
<specular.ode.result.ODEResult at 0x1765982d8d0>

To access the numerical results, call .history(). It returns a tuple containing the time grid and the numerical solution.

import specular

def F(t, u):
    return -2*u 

specular.Euler_scheme(of_Type=1, F=F, t_0=0.0, u_0=1.0, T=2.5, h=0.1).history()

Running the specular Euler scheme of Type 1: 100%|██████████| 24/24 [00:00<?, ?it/s]
(array([0. , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1. , 1.1, 1.2,
        1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2. , 2.1, 2.2, 2.3, 2.4, 2.5]),
 array([1.        , 0.8       , 0.62169432, 0.48101574, 0.37172557,
        0.2870388 , 0.22149069, 0.17081087, 0.13166787, 0.1014624 ,
        0.07816953, 0.06021577, 0.04638162, 0.0357239 , 0.02751427,
        0.02119088, 0.01632056, 0.0125695 , 0.00968054, 0.00745555,
        0.00574195, 0.00442221, 0.00340579, 0.00262299, 0.00202011,
        0.0015558 ]))

To visualize the numerical results, call .visualization().

import specular
import numpy as np

def F(t, u):
   return -2*u 

def exact_sol(t):
    return np.exp(-2*t)

def u_0(t_0):
    return exact_sol(t_0)

specular.Euler_scheme(of_Type='1', F=F, t_0=0.0, u_0=u_0, T=2.5, h=0.1).visualization(exact_sol=exact_sol, save_path="specular-Euler-scheme-of-Type-1")

Running the specular Euler scheme of Type 1: 100%|██████████| 24/24 [00:00<?, ?it/s]
Figure saved: figures\specular-Euler-scheme-of-Type-1

specular-Euler-scheme-of-Type-1

To obtain the table of the numerical results, call .table(). 

import specular
import numpy as np

def F(t, u):
    return -2*u

def exact_sol(t):
    return np.exp(-2*t)

def u_0(t_0):
    return exact_sol(t_0)

specular.Euler_scheme(of_Type=4, F=F, t_0=0.0, u_0=u_0, T=2.5, h=0.1).table(exact_sol=exact_sol, save_path="specular-Euler-scheme-of-type-4")

Running the specular Euler scheme of Type 4: 100%|██████████| 25/25 [00:00<?, ?it/s]
Table saved: tables\specular-Euler-scheme-of-type-4.csv

.visualization() and .table() are chainable.

import specular
import numpy as np

def F(t, u):
    return -2*u

def exact_sol(t):
    return np.exp(-2*t)

def u_0(t_0):
    return exact_sol(t_0)

specular.Euler_scheme(of_Type=4, F=F, t_0=0.0, u_0=u_0, T=2.5, h=0.1).visualization(exact_sol=exact_sol).table(exact_sol=exact_sol)

Running the specular Euler scheme of Type 4: 100%|██████████| 25/25 [00:00<?, ?it/s]

To compute the total error of the numerical results, call .total_error(). The exact solution is required. The norm can be max, l1, or l2.

def F(t, u):
    return -2*u 

def exact_sol(t):
    return np.exp(-2*t)

def u_0(t_0):
    return exact_sol(t_0)

specular.Euler_scheme(of_Type=5, F=F, t_0=0.0, u_0=u_0, T=10.0, h=0.1).total_error(exact_sol=exact_sol, norm='max')

Running the specular Euler scheme of Type 5: 100%|██████████| 100/100 [00:00<00:00, 300882.64it/s]
0.0011409613137273178

2.2.2. Specular trigonometric scheme

import specular

def F(t, u):
    return -2*u 

def exact_sol(t):
    return np.exp(-2*t)

def u_0(t_0):
    return exact_sol(t_0)

t_0 = 0.0
h = 0.1
u_1 = exact_sol(t_0 + h)

specular.trigonometric_scheme(F=F, t_0=t_0, u_0=u_0, u_1=u_1, T=2.5, h=h).visualization(exact_sol=exact_sol, save_path="specular-trigonometric")

Running specular trigonometric scheme: 100%|██████████| 24/24 [00:00<?, ?it/s]
Figure saved: figures\specular-trigonometric

specular-trigonometric-scheme

2.2.3. Classical schemes

The three classical schemes are available: the explicit Euler, the implicit Euler, and the Crank-Nicolson schemes.

import specular
import numpy as np
import matplotlib.pyplot as plt

def F(t, u):
    return -(t*u)/(1-t**2)
def exact_sol(t):
    return np.sqrt(1 - t**2)
def u_0(t_0):
    return exact_sol(t_0)
t_0 = 0.0
T = 0.9
h = 0.05

result_EE = specular.classical_scheme(F=F, t_0=t_0, u_0=u_0, T=T, h=h, form="explicit Euler").history()
result_IE = specular.classical_scheme(F=F, t_0=t_0, u_0=u_0, T=T, h=h, form="implicit Euler").history()
result_CN = specular.classical_scheme(F=F, t_0=t_0, u_0=u_0, T=T, h=h, form="Crank-Nicolson").history()
exact_values = np.array([exact_sol(t) for t in result_EE[0]])

plt.figure(figsize=(5.5, 2.5))

plt.plot(result_EE[0], exact_values, color='black', label='Exact solution')
plt.plot(result_EE[0], result_EE[1],  marker='x', linestyle='None', markerfacecolor='none', markeredgecolor='red', label='Explicit Euler') 
plt.plot(result_IE[0], result_IE[1],  marker='x', linestyle='None', markerfacecolor='none', markeredgecolor='blue', label='Implicit Euler') 
plt.plot(result_CN[0], result_CN[1],  marker='x', linestyle='None', markerfacecolor='none', markeredgecolor='purple', label='Crank-Nicolson')

plt.xlabel(r"Time", fontsize=10)
plt.ylabel(r"Solution", fontsize=10)
plt.grid(True)
plt.legend(loc='center left', bbox_to_anchor=(1.02, 0.5), borderaxespad=0., fontsize=10)
plt.savefig('figures/classical-schemes.png', dpi=1000, bbox_inches='tight')
plt.show()

Running the explicit Euler scheme: 100%|██████████| 18/18 [00:00<?, ?it/s]
Running the implicit Euler scheme: 100%|██████████| 18/18 [00:00<?, ?it/s]
Running Crank-Nicolson scheme: 100%|██████████| 18/18 [00:00<00:00, 17988.44it/s]

classical-schemes

2.2.4. API Reference

specular.ode.solver

Let the source function \(F:[t_0, T] \times \mathbb{R} \to \mathbb{R}\) be given, and the initial data \(u_0:\mathbb{R} \to \mathbb{R}\) be given. Consider the initial value problem:

\[ u'(t) = F(t, u(t)) \qquad \text{(IVP)} \]

with the initial condition \(u(t_0) = u_0(t_0)\). To solve (IVP) numerically, this module provides implementations of the specular Euler schemes and the specular trigonometric scheme.

Euler_scheme(of_Type, F, t_0, u_0, T, h=1e-06, u_1=None, tol=1e-12, zero_tol=1e-08, max_iter=100)

Solves an initial value problem (IVP) using the specular Euler scheme of Type 1, 2, 3, 4, 5, and 6.

Parameters:

Name Type Description Default
of_Type int | str

The type of the specular Euler scheme. Options: 1, '1', 2, '2', 3, '3', 4, '4', 5, '5', 6, '6'.

 required
F callable

The given source function F in (IVP). The calling signature should be F(t, u) where t and u are scalars.

 required
t_0 float

The starting time of the simulation.

 required
u_0 callable

The given initial condition u_0 in (IVP).

 required
T float

The end time of the simulation.

 required
h float

Mesh size used in the finite difference approximation. Must be positive.

 1e-06
u_1 callable | float | None

The numerical solution at the time t_1 = t_0 + h. If a float or callable is provided, it is used as the first-step value.

 None
tol float | optional

Tolerance for fixed-point iteration Used for Types 3, 4, 5, and 6.

 1e-12
zero_tol float | floating

A small threshold used to determine if the denominator (alpha + beta) is close to zero for numerical stability.

 1e-08
max_iter int | optional

Max iterations for fixed-point solver.

 100

Returns:

Type Description
ODEResult

An object containing (t, u) data and the scheme name.
Source code in specular\ode\solver.py 
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
def Euler_scheme(
    of_Type: int | str,
    F: Callable[[float, float], float],
    t_0: float, 
    u_0: Callable[[float], float] | float,
    T: float, 
    h: float = 1e-6,
    u_1: Callable[[float], float] | float | None = None,
    tol: float = 1e-12, 
    zero_tol: float = 1e-8,
    max_iter: int = 100
) -> ODEResult:
    """
    Solves an initial value problem (IVP) using the specular Euler scheme of Type 1, 2, 3, 4, 5, and 6.

    Parameters:
        of_Type (int | str):
            The type of the specular Euler scheme.
            Options: ``1``, ``'1'``, ``2``, ``'2'``, ``3``, ``'3'``, ``4``, ``'4'``, ``5``, ``'5'``, ``6``, ``'6'``.
        F (callable):
            The given source function ``F`` in (IVP).
            The calling signature should be ``F(t, u)`` where ``t`` and ``u`` are scalars.
        t_0 (float):
            The starting time of the simulation.
        u_0 (callable):
            The given initial condition ``u_0`` in (IVP).
        T (float):
            The end time of the simulation.
        h (float, optional):
            Mesh size used in the finite difference approximation. Must be positive.
        u_1 (callable | float | None):
            The numerical solution at the time ``t_1 = t_0 + h``.
            If a float or callable is provided, it is used as the first-step value.
        tol (float | optional):
            Tolerance for fixed-point iteration
            Used for Types 3, 4, 5, and 6.
        zero_tol (float | np.floating):
            A small threshold used to determine if the denominator (alpha + beta) is close to zero for numerical stability.
        max_iter (int | optional):
            Max iterations for fixed-point solver.

    Returns:
        An object containing ``(t, u)`` data and the scheme name.
    """
    Type = str(of_Type)

    scheme = 'specular Euler scheme of Type ' + Type
    steps = _num_steps(t_0, T, h)

    all_history = {}
    t_history = []
    u_history = []

    if Type in ['1', '2', '3']:
        t_prev = t_0
        u_prev = u_0(t_0) if callable(u_0) else u_0

        t_history.append(t_prev)
        u_history.append(u_prev)

        t_curr = t_prev + h

        if u_1 is None:
            # explicit Euler to get u_1
            u_curr = u_prev + h * F(t_prev, u_prev)
        elif callable(u_1):
            u_curr = u_1(t_curr)
        else:
            u_curr = u_1

        t_history.append(t_curr)
        u_history.append(u_curr)

        if Type == '1':
            _of_Type_1(F, h, zero_tol, steps, t_prev, t_curr, u_prev, u_curr, t_history, u_history)

        elif Type == '2':
            _of_Type_2(F, h, zero_tol, steps, t_prev, t_curr, u_prev, u_curr, t_history, u_history)

        elif Type == '3':
            _of_Type_3(F, h, tol, zero_tol, max_iter, steps, t_prev, t_curr, u_prev, u_curr, t_history, u_history)

    elif Type in ['4', '5', '6']:
        t_curr = t_0
        u_curr = u_0(t_0) if callable(u_0) else u_0  

        t_history.append(t_curr)
        u_history.append(u_curr)

        if Type == '4':
            _of_Type_4(F, h, tol, zero_tol, max_iter, steps, t_curr, u_curr, t_history, u_history)

        elif Type == '5':
            _of_Type_5(F, h, tol, zero_tol, max_iter, steps, t_curr, u_curr, t_history, u_history)

        elif Type == '6':
            _of_Type_6(F, h, tol, zero_tol, max_iter, steps, t_curr, u_curr, t_history, u_history)

    else:
        raise ValueError(f"Unknown type. Got {of_Type}. Supported types: '1', '2', '3', '4', '5', and '6'")

    all_history["variables"] = np.array(t_history)
    all_history["values"] = np.array(u_history)

    return ODEResult(
        scheme=scheme,
        h=h,
        all_history=all_history
    )

Heun_scheme(F, t_0, u_0, T, h=1e-06, zero_tol=1e-08)

Solves an initial value problem using the specular Heun scheme.

This is the one-step, two-stage scheme obtained from one fixed-point iteration of the specular Euler scheme of Type 6.

Source code in specular\ode\solver.py 
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
def Heun_scheme(
    F: Callable[[float, float], float],
    t_0: float,
    u_0: Callable[[float], float] | float,
    T: float,
    h: float = 1e-6,
    zero_tol: float = 1e-8,
) -> ODEResult:
    """
    Solves an initial value problem using the specular Heun scheme.

    This is the one-step, two-stage scheme obtained from one fixed-point iteration of the specular Euler scheme of Type 6.
    """
    t_curr = t_0
    u_curr = u_0(t_0) if callable(u_0) else u_0

    all_history = {}
    t_history = [t_curr]
    u_history = [u_curr]

    steps = _num_steps(t_0, T, h)

    for _ in tqdm(range(steps), desc="Running the specular Heun scheme"):
        t_next = t_curr + h

        beta = F(t_curr, u_curr)
        u_predictor = u_curr + h * beta
        alpha = F(t_next, u_predictor)

        u_next = u_curr + h * A(alpha, beta, zero_tol=zero_tol)

        t_curr, u_curr = t_next, u_next

        t_history.append(t_curr)
        u_history.append(u_curr)

    all_history["variables"] = np.array(t_history)
    all_history["values"] = np.array(u_history)

    return ODEResult(
        scheme="specular Heun scheme",
        h=h,
        all_history=all_history,
    )

ellipse_scheme(F, t_0, u_0, T, a, b, h=1e-06, tol=1e-12, zero_tol=1e-08, max_iter=100)

Solves an initial value problem (IVP) using the specular ellipse scheme.

Parameters:

Name Type Description Default
F callable

The given source function F in (IVP). The calling signature should be F(t, u) where t and u are scalars.

 required
t_0 float

The starting time of the simulation.

 required
u_0 callable | float

The given initial condition u_0 in (IVP). If a callable is provided, the initial value is determined by calling u_0(t_0).

 required
T float

The end time of the simulation.

 required
a float

The semi-axis length in the t direction. Must be positive.

 required
b float

The semi-axis length in the u direction. Must be positive.

 required
h float

Mesh size used in the finite difference approximation. Must be positive.

 1e-06
tol float

Tolerance for fixed-point iteration.

 1e-12
zero_tol float | floating

A small threshold used to determine if the denominator (alpha + beta) is close to zero for numerical stability.

 1e-08
max_iter int

Max iterations for fixed-point solver.

 100

Returns:

Type Description
ODEResult

An object containing (t, u) data and the scheme name.
Source code in specular\ode\solver.py 
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
def ellipse_scheme(
    F: Callable[[float, float], float],
    t_0: float,
    u_0: Callable[[float], float] | float,
    T: float,
    a: float,
    b: float,
    h: float = 1e-6,
    tol: float = 1e-12,
    zero_tol: float = 1e-8,
    max_iter: int = 100,
) -> ODEResult:
    """
    Solves an initial value problem (IVP) using the specular ellipse scheme.

    Parameters:
        F (callable):
            The given source function ``F`` in (IVP).
            The calling signature should be ``F(t, u)`` where ``t`` and ``u`` are scalars.
        t_0 (float):
            The starting time of the simulation.
        u_0 (callable | float):
            The given initial condition ``u_0`` in (IVP).
            If a callable is provided, the initial value is determined by calling ``u_0(t_0)``.
        T (float):
            The end time of the simulation.
        a (float):
            The semi-axis length in the ``t`` direction. Must be positive.
        b (float):
            The semi-axis length in the ``u`` direction. Must be positive.
        h (float, optional):
            Mesh size used in the finite difference approximation. Must be positive.
        tol (float, optional):
            Tolerance for fixed-point iteration.
        zero_tol (float | np.floating):
            A small threshold used to determine if the denominator (alpha + beta) is close to zero for numerical stability.
        max_iter (int, optional):
            Max iterations for fixed-point solver.

    Returns:
        An object containing ``(t, u)`` data and the scheme name.
    """
    if a <= 0 or b <= 0:
        raise ValueError(f"Semi-axes 'a' and 'b' must be positive. Got a={a}, b={b}")

    t_curr = t_0
    u_curr = u_0(t_0) if callable(u_0) else u_0

    all_history = {}
    t_history = [t_curr]
    u_history = [u_curr]

    steps = _num_steps(t_0, T, h)

    scale = b / a
    inv_scale = a / b

    for k in tqdm(range(steps), desc="Running the specular ellipse scheme"):
        beta = inv_scale * F(t_curr, u_curr)
        t_next = t_curr + h

        # Initial guess: explicit Euler
        u_temp = u_curr + h * F(t_curr, u_curr)
        u_guess = u_temp

        # Fixed-point iteration
        for _ in range(max_iter):
            alpha = inv_scale * F(t_next, u_temp)
            A_val = float(A(alpha, beta, zero_tol=zero_tol))
            u_guess = u_curr + h * scale * A_val

            if abs(u_guess - u_temp) < tol:
                break

            u_temp = u_guess
        else:
            print(f"Warning: fixed-point iteration did not converge at step {k+1}")

        t_curr, u_curr = t_next, u_guess

        t_history.append(t_curr)
        u_history.append(u_curr)

    all_history["variables"] = np.array(t_history)
    all_history["values"] = np.array(u_history)

    return ODEResult(
        scheme="specular ellipse scheme",
        h=h,
        all_history=all_history,
    )

trigonometric_scheme(F, t_0, u_0, u_1, T, h=1e-06)

Solves an initial value problem (IVP) using the specular trigonometric scheme.

Parameters:

Name Type Description Default
F callable

The given source function F in (IVP). The calling signature should be F(t, u) where t and u are scalars.

 required
t_0 float

The starting time of the simulation.

 required
u_0 callable

The given initial condition u_0 in (IVP).

 required
u_1 callable | float

The numerical solution at the time t_1 = t_0 + h. If a float or callable is provided, it is used as the first-step value.

 required
T float

The end time of the simulation.

 required
h float

Mesh size used in the finite difference approximation. Must be positive.

 1e-06

Returns:

Type Description
ODEResult

An object containing (t, u) data and the scheme name.
Source code in specular\ode\solver.py 
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
def trigonometric_scheme(
    F: Callable[[float, float], float],
    t_0: float, 
    u_0: Callable[[float], float] | float,
    u_1: Callable[[float], float] | float,
    T: float, 
    h: float = 1e-6
) -> ODEResult:
    """
    Solves an initial value problem (IVP) using the specular trigonometric scheme.

    Parameters:
        F (callable):
            The given source function ``F`` in (IVP).
            The calling signature should be ``F(t, u)`` where ``t`` and ``u`` are scalars.
        t_0 (float):
            The starting time of the simulation.
        u_0 (callable):
            The given initial condition ``u_0`` in (IVP).
        u_1 (callable | float):
            The numerical solution at the time ``t_1 = t_0 + h``.
            If a float or callable is provided, it is used as the first-step value.
        T (float):
            The end time of the simulation.
        h (float, optional):
            Mesh size used in the finite difference approximation. Must be positive.

    Returns:
        An object containing ``(t, u)`` data and the scheme name.
    """
    t_prev = t_0
    u_prev = u_0(t_0) if callable(u_0) else u_0

    t_curr = t_0 + h
    u_curr = u_1(t_curr) if callable(u_1) else u_1

    all_history = {}
    t_history = [t_prev, t_curr]
    u_history = [u_prev, u_curr]

    steps = _num_steps(t_0, T, h)

    for _ in tqdm(range(steps - 1), desc="Running specular trigonometric scheme"):
        t_next = t_curr + h
        u_next = u_curr + h * math.tan(2 * math.atan(F(t_curr, u_curr)) - math.atan((u_curr - u_prev) / h))

        t_history.append(t_next)
        u_history.append(u_next)

        t_prev, u_prev = t_curr, u_curr
        t_curr, u_curr = t_next, u_next

    all_history["variables"] = np.array(t_history)
    all_history["values"] = np.array(u_history)

    return ODEResult(
        scheme="specular trigonometric scheme",
        h=h,
        all_history=all_history
    )

specular.ode.classical_solver

Classical fixed-step schemes for scalar first-order ODEs.

specular.ode.result

ODEResult

Source code in specular\ode\result.py 
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
class ODEResult:
    def __init__(
        self,
        scheme: str,
        h: float,
        all_history: dict
    ):
        self.scheme = scheme
        self.h = h
        self.time_grid = all_history["variables"]
        self.numerical_sol = all_history["values"]

    def history(
        self
    ) -> Tuple[np.ndarray, np.ndarray]:
        """
        Returns the time grid and the numerical solution as a tuple.

        Returns:
            (time_grid, numerical_sol)
        """
        return self.time_grid, self.numerical_sol

    def visualization(
        self, 
        figure_size: tuple = (5.5, 2.5), 
        exact_sol: Optional[Callable[[float], float]] = None,
        save_path: Optional[str] = None
    ):
        import matplotlib.pyplot as plt

        plt.figure(figsize=figure_size)

        if exact_sol is not None:
            exact_values = np.array([exact_sol(t) for t in self.time_grid])
            plt.plot(self.time_grid, exact_values, color='black', label='Exact solution')

        number_of_circles = max(1, len(self.time_grid) // 30)

        capitalized_name = self.scheme[0].upper() + self.scheme[1:]

        plt.plot(self.time_grid, self.numerical_sol, linestyle='--', marker='o', color='red', markersize=5, markevery=number_of_circles, markerfacecolor='none', markeredgewidth=1.0, label=capitalized_name)

        plt.xlabel(r"Time", fontsize=10)
        plt.ylabel(r"Solution", fontsize=10)
        plt.grid(True)
        plt.legend(loc='center left', bbox_to_anchor=(1.02, 0.5), borderaxespad=0., fontsize=10)

        if save_path:
            if not os.path.exists('figures'):
                os.makedirs('figures')

            save_path = save_path.replace(" ", "-")
            full_path = os.path.join("figures", save_path)

            if not save_path.endswith(".png"):
                save_path += ".png"

            plt.savefig(full_path, dpi=1000, bbox_inches='tight')

            print(f"Figure saved: {full_path}")

        plt.show()

        return self

    def table(self,
        exact_sol: Optional[Callable[[float], float]] = None,
        save_path: Optional[str] = None
    ):
        import pandas as pd

        result = pd.DataFrame(self.numerical_sol, index=self.time_grid, columns=["Numerical solution"])
        result.index.name = "Time"

        if exact_sol:
            result["Exact solution"] = [exact_sol(t) for t in self.time_grid]
            result["Error"] = abs(result["Numerical solution"] - result["Exact solution"])

        if save_path:
            if not os.path.exists('tables'):
                os.makedirs('tables')

            save_path = save_path.replace(" ", "-")
            full_path = os.path.join("tables", save_path)

            if full_path.endswith(".txt"):
                with open(full_path, "w") as f:
                    f.write(result.to_string())
            else:
                if not full_path.endswith(".csv"):
                    full_path += ".csv"

                result.to_csv(full_path)

            print(f"Table saved: {full_path}")

        return result

    def total_error(self,
        exact_sol: Callable[[float], float] | list | np.ndarray,
        norm: str = 'max'
    ) -> float:
        """
        Calculates the error between the numerical solution and the exact solution.

        Parameters:
            exact_sol (callable | list | np.ndarray):
                A function that returns the exact solution at a given time ``t``, or a list/array containing the exact values corresponding to ``time_grid``.
            norm (str | optional):
                The type of norm to use ``'max'``, ``'l2'``, or ``'l1'``.

        Returns:
            The computed error value.

        Raises:
            TypeError:
                If ``exact_sol`` is neither a callable nor a list/array.
            ValueError:
                If ``exact_sol`` (list) shape does not match ``numerical_sol``.
        """
        if callable(exact_sol):
            exact_values = np.array([exact_sol(t) for t in self.time_grid])

        elif isinstance(exact_sol, (list, np.ndarray)):
            exact_values = np.asarray(exact_sol, dtype=float) 

        else:
            raise TypeError("exact_sol must be a callable or a list/array.")

        if exact_values.shape != self.numerical_sol.shape:
             raise ValueError(f"Shape mismatch: exact_sol {exact_values.shape} vs numerical_sol {self.numerical_sol.shape}")

        error_vector = np.abs(exact_values - self.numerical_sol)

        if norm == 'max':
            return float(np.max(error_vector))

        elif norm == 'l2':
            return float(np.sqrt(np.sum(error_vector**2) * self.h))

        elif norm == 'l1':
            return float(np.sum(error_vector) * self.h)

        else:
            raise ValueError(f"Unknown norm type. Got '{norm}'. Supported types: 'max', 'l2', 'l1'.")

history()

Returns the time grid and the numerical solution as a tuple.

Returns:

Type Description
Tuple[ndarray, ndarray]

(time_grid, numerical_sol)
Source code in specular\ode\result.py 
17
18
19
20
21
22
23
24
25
26
def history(
    self
) -> Tuple[np.ndarray, np.ndarray]:
    """
    Returns the time grid and the numerical solution as a tuple.

    Returns:
        (time_grid, numerical_sol)
    """
    return self.time_grid, self.numerical_sol

total_error(exact_sol, norm='max')

Calculates the error between the numerical solution and the exact solution.

Parameters:

Name Type Description Default
exact_sol callable | list | ndarray

A function that returns the exact solution at a given time t, or a list/array containing the exact values corresponding to time_grid.

 required
norm str | optional

The type of norm to use 'max', 'l2', or 'l1'.

 'max'

Returns:

Type Description
float

The computed error value.

Raises:

Type Description
TypeError

If exact_sol is neither a callable nor a list/array.
ValueError

If exact_sol (list) shape does not match numerical_sol.
Source code in specular\ode\result.py 
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
def total_error(self,
    exact_sol: Callable[[float], float] | list | np.ndarray,
    norm: str = 'max'
) -> float:
    """
    Calculates the error between the numerical solution and the exact solution.

    Parameters:
        exact_sol (callable | list | np.ndarray):
            A function that returns the exact solution at a given time ``t``, or a list/array containing the exact values corresponding to ``time_grid``.
        norm (str | optional):
            The type of norm to use ``'max'``, ``'l2'``, or ``'l1'``.

    Returns:
        The computed error value.

    Raises:
        TypeError:
            If ``exact_sol`` is neither a callable nor a list/array.
        ValueError:
            If ``exact_sol`` (list) shape does not match ``numerical_sol``.
    """
    if callable(exact_sol):
        exact_values = np.array([exact_sol(t) for t in self.time_grid])

    elif isinstance(exact_sol, (list, np.ndarray)):
        exact_values = np.asarray(exact_sol, dtype=float) 

    else:
        raise TypeError("exact_sol must be a callable or a list/array.")

    if exact_values.shape != self.numerical_sol.shape:
         raise ValueError(f"Shape mismatch: exact_sol {exact_values.shape} vs numerical_sol {self.numerical_sol.shape}")

    error_vector = np.abs(exact_values - self.numerical_sol)

    if norm == 'max':
        return float(np.max(error_vector))

    elif norm == 'l2':
        return float(np.sqrt(np.sum(error_vector**2) * self.h))

    elif norm == 'l1':
        return float(np.sum(error_vector) * self.h)

    else:
        raise ValueError(f"Unknown norm type. Got '{norm}'. Supported types: 'max', 'l2', 'l1'.")