Source code for simulai.utilities.problem_classes

# (C) Copyright IBM Corp. 2019, 2020, 2021, 2022.

#    Licensed under the Apache License, Version 2.0 (the "License");
#    you may not use this file except in compliance with the License.
#    You may obtain a copy of the License at

#           http://www.apache.org/licenses/LICENSE-2.0

#     Unless required by applicable law or agreed to in writing, software
#     distributed under the License is distributed on an "AS IS" BASIS,
#     WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#     See the License for the specific language governing permissions and
#     limitations under the License.

import numpy as np
from scipy import linalg


[docs]class NonlinearOscillator: def __init__(self, forcing : bool =False, p : int =3, alpha1=-0.1, alpha2=-2, beta1=2, beta2=-0.1): self.alpha1 = alpha1 self.alpha2 = alpha2 self.beta1 = beta1 self.beta2 = beta2 self.p = p if forcing is True: self.eval = self._eval_with_forcing elif forcing is False: self.eval = self._eval_without_forcing else: raise Exception(f"It is expected forcing to be bool but received {type(forcing)}") def _eval_without_forcing(self, state): state_ = state[0] x = state_[0] y = state_[1] f = self.alpha1 * (x ** self.p) + self.beta1 * (y ** self.p) g = self.alpha2 * (x ** self.p) + self.beta2 * (y ** self.p) return f, g def _eval_with_forcing(self, state): state_ = state[0] x = state_[0] y = state_[1] u = state_[2] v = state_[3] f = self.alpha1 * (x ** self.p) + self.beta1 * (y ** self.p) + u g = self.alpha2 * (x ** self.p) + self.beta2 * (y ** self.p) + v return f, g def __call__(self, state): f, g = self.eval(state) return np.array([f, g])
[docs]class LorenzSystem: def __init__(self, rho, sigma, beta, forcing : bool =False, use_t : bool =False): """ :param rho: :param sigma: :param beta: """ self.rho = rho self.beta = beta self.sigma = sigma if forcing is True: self.eval = self._eval_with_forcing elif forcing is False and use_t is False: self.eval = self._eval_without_forcing elif forcing is False and use_t is True: self.eval = self._eval_without_forcing_t else: raise Exception(f"It is expected forcing to be bool but received {type(forcing)}") def _eval_without_forcing(self, state, *args): state_ = state[0] x = state_[0] y = state_[1] z = state_[2] f = self.sigma * (y - x) g = x*(self.rho - z) - y h = x*y - self.beta*z return np.array([f, g, h]) def _eval_without_forcing_t(self, t, state): state_ = state x = state_[0] y = state_[1] z = state_[2] f = self.sigma * (y - x) g = x*(self.rho - z) - y h = x*y - self.beta*z return np.array([f, g, h]) def _eval_with_forcing(self, state): state_ = state[0] x = state_[0] y = state_[1] z = state_[2] u = state_[3] v = state_[4] w = state_[5] f = self.sigma * (y - x) + u g = x*(self.rho - z) - y + v h = x*y - self.beta*z + w return np.array([f, g, h]) def __call__(self, state): f, g, h = self.eval(state) return np.array([f, g, h])
[docs] def jacobian(self, state, e, w, dt): x = state[0] y = state[1] z = state[2] e1 = e[0] e2 = e[1] e3 = e[2] D = np.array([ [-self.sigma, self.sigma, 0], [-z + self.rho, -1, -x], [y, x, -self.beta] ]) J = np.eye(3) + dt*D w_prev = w w = linalg.orth(J*w_prev) de1 = np.log(np.linalg.norm(w[:, 0], 2)) de2 = np.log(np.linalg.norm(w[:, 1], 2)) de3 = np.log(np.linalg.norm(w[:, 2], 2)) e1 += de1 e2 += de2 e3 += de3 w1 = w[:, 0] / np.linalg.norm(w[:, 0], 2) w2 = w[:, 1] / np.linalg.norm(w[:, 1], 2) w3 = w[:, 2] / np.linalg.norm(w[:, 2], 2) return np.array([e1, e2, e3]), np.hstack([w1[:, None], w2[:, None], w3[:, None]])