Coverage for pySDC/projects/parallelSDC_reloaded/lorenz_setup.py: 100%

1#!/usr/bin/env python3 

2# -*- coding: utf-8 -*- 

3""" 

4Created on Tue Dec 5 10:08:04 2023 

5 

6Script to numerically compute number revolution periods for the Lorenz system 

7""" 

8import numpy as np 

9from scipy import signal 

10import matplotlib.pyplot as plt 

11 

12from pySDC.projects.parallelSDC_reloaded.utils import solutionExact 

13 

14script = __file__.split('/')[-1].split('.')[0] 

15 

16tEnd = 2 

17nSteps = tEnd * 50 

18tVals = np.linspace(0, tEnd, nSteps + 1) 

19 

20nPeriods = 2 

21 

22print(f"Computing exact solution up to t={tEnd} ...") 

23uExact = solutionExact(tEnd, nSteps, "LORENZ", u0=(5, -5, 20)) 

24 

25z = uExact[:, -1] 

26idx = signal.find_peaks(z)[0][nPeriods - 1] 

27 

28 

29print(f'tEnd for {nPeriods} periods : {tVals[idx]}') 

30 

31figName = f"{script}_traj" 

32 

33plt.figure(figName) 

34plt.plot(tVals, uExact[:, 0], '-', label="$x(t)$") 

35plt.plot(tVals, uExact[:, 1], '-', label="$y(t)$") 

36plt.plot(tVals, uExact[:, 2], '-', label="$z(t)$") 

37plt.vlines(tVals[idx], ymin=-20, ymax=40, linestyles="--", linewidth=1) 

38 

39plt.legend() 

40plt.xlabel("time") 

41plt.ylabel("trajectory") 

42plt.tight_layout()