#!/usr/bin/env python # -*- coding: utf-8 -*- """ *********************************************************************************** tutorial18.py DAE Tools: pyDAE module, www.daetools.com Copyright (C) Dragan Nikolic *********************************************************************************** DAE Tools is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License version 3 as published by the Free Software Foundation. DAE Tools is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with the DAE Tools software; if not, see <http://www.gnu.org/licenses/>. ************************************************************************************ """ __doc__ = """ This tutorial shows one more problem solved using the NumPy arrays that operate on DAE Tools variables. The model is taken from the Sundials ARKODE (ark_analytic_sys.cpp). The ODE system is defined by the following system of equations: .. code-block:: none dy/dt = A*y where: .. code-block:: none A = V * D * Vi V = [1 -1 1; -1 2 1; 0 -1 2]; Vi = 0.25 * [5 1 -3; 2 2 -2; 1 1 1]; D = [-0.5 0 0; 0 -0.1 0; 0 0 lam]; lam is a large negative number. The analytical solution to this problem is: .. code-block:: none Y(t) = V*exp(D*t)*Vi*Y0 for t in the interval [0.0, 0.05], with initial condition y(0) = [1,1,1]'. The stiffness of the problem is directly proportional to the value of "lamda". The value of lamda should be negative to result in a well-posed ODE; for values with magnitude larger than 100 the problem becomes quite stiff. In this example, we choose lamda = -100. The solution: .. code-block:: none lamda = -100 reltol = 1e-06 abstol = 1e-10 -------------------------------------- t y0 y1 y2 -------------------------------------- 0.0050 0.70327 0.70627 0.41004 0.0100 0.52267 0.52865 0.05231 0.0150 0.41249 0.42145 -0.16456 0.0200 0.34504 0.35696 -0.29600 0.0250 0.30349 0.31838 -0.37563 0.0300 0.27767 0.29551 -0.42383 0.0350 0.26138 0.28216 -0.45296 0.0400 0.25088 0.27459 -0.47053 0.0450 0.24389 0.27053 -0.48109 0.0500 0.23903 0.26858 -0.48740 -------------------------------------- The plot of the 'y0', 'y1', 'y2' variables: .. image:: _static/tutorial18-results.png :width: 500px """ import sys, numpy, scipy.linalg from time import localtime, strftime from daetools.pyDAE import * # Standard variable types are defined in variable_types.py from pyUnits import m, kg, s, K, Pa, mol, J, W typeNone = daeVariableType("typeNone", unit(), 0, 1E10, 0, 1e-10) class modTutorial(daeModel): def __init__(self, Name, Parent = None, Description = ""): daeModel.__init__(self, Name, Parent, Description) self.x = daeDomain("x", self, unit(), "") self.y = daeVariable("y", typeNone, self, "", [self.x]) def DeclareEquations(self): daeModel.DeclareEquations(self) # Input data: lamda = -100 V = numpy.array([[1, -1, 1], [-1, 2, 1], [0, -1, 2]]) Vi = 0.25 * numpy.array([[5, 1, -3], [2, 2, -2], [1, 1, 1]]) D = numpy.array([[-0.5, 0, 0], [0, -0.1, 0], [0, 0, lamda]]) self.y0 = numpy.array([1.0, 1.0, 1.0]) # Create a vector of y: y = numpy.empty(3, dtype=object) y[:] = [self.y(i) for i in range(3)] # Create a vector of dy/dt: dydt = numpy.empty(3, dtype=object) dydt[:] = [dt(self.y(i)) for i in range(3)] # Create the ODE system: dy/dt = A*y # Use dot product (numpy arrays don't behave as matrices) # or use numpy.matrix where the operator * performs the dot product. Ay = V.dot(D).dot(Vi).dot(y) #print(Ay) for i in range(3): eq = self.CreateEquation("y(%d)" % i, "") eq.Residual = dydt[i] - Ay[i] eq.CheckUnitsConsistency = False class simTutorial(daeSimulation): def __init__(self): daeSimulation.__init__(self) self.m = modTutorial("tutorial18") self.m.Description = __doc__ def SetUpParametersAndDomains(self): self.m.x.CreateArray(3) def SetUpVariables(self): self.m.y.SetInitialConditions(self.m.y0) def run(**kwargs): simulation = simTutorial() relativeTolerance = 1E-6 return daeActivity.simulate(simulation, reportingInterval = 0.005, timeHorizon = 0.050, relativeTolerance = relativeTolerance, **kwargs) if __name__ == "__main__": guiRun = False if (len(sys.argv) > 1 and sys.argv[1] == 'console') else True run(guiRun = guiRun)