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

8 

9import numpy as np 

10from scipy import signal 

11import matplotlib.pyplot as plt 

12 

13from pySDC.projects.parallelSDC_reloaded.utils import solutionExact 

14 

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

16 

17tEnd = 2 

18nSteps = tEnd * 50 

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

20 

21nPeriods = 2 

22 

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

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

25 

26z = uExact[:, -1] 

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

28 

29 

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

31 

32figName = f"{script}_traj" 

33 

34plt.figure(figName) 

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

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

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

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

39 

40plt.legend() 

41plt.xlabel("time") 

42plt.ylabel("trajectory") 

43plt.tight_layout()