Hi guys, I'm working on a larger project, solving non-linear differential equation systems by calculating the steady states. I want to create a very flexible script were I enter a n-dimensional equation system (mostly n=3). The script calculates the steady states and does a stabilty analysis and will have a latex file output.
The main part of the script is working. (see test.py or test.sh for compiling the latex code) But I face two problems I can't solve. My first problem is also discussed here: https://groups.google.com/forum/#!topic/sympy/ebGJLKeMODQ All Variables of the System (here X,Y,Z) should be >= 0. If you look at the last steady state in the latex file, the expressions for Y,Z are quite uncomfortable because of the negative sign. (1, -(\RXZ*\RZY + \RZY + 1)/(\RYZ*\RZY - 1), -(\RXZ + \RYZ + 1)/(\RYZ*\RZY - 1)) So I would like to replace these expression by (1, (\RXZ*\RZY + \RZY + 1)/(1 - \RYZ*\RZY ), (\RXZ + \RYZ + 1)/(1 - \RYZ*\RZY )) If you also have a look at the diagonal elements J_2,2 and J_3,3 of the jacobian matrix. These elements could be replaced by a much shorter term like -\RY*\Y for J_2,2 and -\RZ*\Z for J_3,3. But the function sympy.subs() doesn't work in this case because of the negative sign of the solutions of Y and Z I think I could solve this problem if I can transform (1, -(\RXZ*\RZY + \RZY + 1)/(\RYZ*\RZY - 1), -(\RXZ + \RYZ + 1)/(\RYZ*\RZY - 1)) to (1, (\RXZ*\RZY + \RZY + 1)/(1 - \RYZ*\RZY ), (\RXZ + \RYZ + 1)/(1 - \RYZ*\RZY )) as written above. Do you know how to transform these lines ? My second problem is the analysis of the calculated eigenvalues of the jacobian matrix. If all eigenvalues are negative, there is a stable steady state. The problem is that the eigenvalues can be a number or an expression. For numbers I can easily ask with an if statement, if they are negative. For the expressions I found the function solve_univariate_inequality() but it doesn't seem very robust if I have a longer expressions with a lot of unknown positive constants. Is there any function which solves ineualities with a flexible amount of input data? Sorry for this long text and thank you in advance for your answers! Alex -- You received this message because you are subscribed to the Google Groups "sympy" group. To unsubscribe from this group and stop receiving emails from it, send an email to [email protected]. To post to this group, send email to [email protected]. Visit this group at https://groups.google.com/group/sympy. To view this discussion on the web visit https://groups.google.com/d/msgid/sympy/59e7ec78-4cda-4ad6-a20c-dc7854b27e67%40googlegroups.com. For more options, visit https://groups.google.com/d/optout.
#!/home/atille/.anaconda2/bin/python #----------------------------------------------------------------------------------------------------------------------# # # # template for calculating steady states of ordinary DES # # and printing the results into a latex file # # # #----------------------------------------------------------------------------------------------------------------------# #----------------------------------------------------------------------------------------------------------------------# # # # required libraries, sympy 1.0 and mpmath # # (on Ubuntu 12.04LTS) sympy 1.0 can be installed in an anaconda2 enviroment # # please look at https://docs.continuum.io/anaconda/install for further explanations # #----------------------------------------------------------------------------------------------------------------------# #----------------------------------------------------------------------------------------------------------------------# # # # the standard latex document requires the .tex file teildok.tex and a preambel # # # #----------------------------------------------------------------------------------------------------------------------# from __future__ import division import sys import mpmath from sympy import * from sympy import solve from sympy import init_printing #init_printing() def print_latex_header(savedata,variables_orig,constants_orig,variables_dimless,constants_dimless,document_type="teildok"): if (document_type=="teildok"): savedata.write(r'\input{teildok}' + '\n') savedata.write(r'\DocumentBegin{\string~/sysbio/arbeit/preambel}' + '\n') else: savedata.write(r'\documentclass[11pt,onehalfspacing,numbers=noenddot,openany]{scrbook}' + '\n') savedata.write(r'\usepackage[utf8]{inputenc}' + '\n') savedata.write(r'\usepackage[utf8]{inputenc}' + '\n') savedata.write(r'\usepackage{oldgerm}' + '\n') savedata.write(r'\usepackage{amsmath,amsthm,amssymb,eurosym,amsfonts}' + '\n') savedata.write(r'\numberwithin{equation}{section}' + '\n') savedata.write(r'\usepackage[ngerman]{babel}' + '\n') savedata.write(r'\usepackage{xspace}' + '\n') savedata.write(r'\usepackage[margin=2cm]{geometry}' + '\n') savedata.write(r'\usepackage{graphicx}' + '\n') savedata.write(r'\usepackage{pdfpages}' + '\n') savedata.write(r'\begin{document}' + '\n') for item in variables_orig: savedata.write(r'\def%s{\mathrm{%s}}' %(latex(item),latex(item)[-1])) savedata.write('\n') savedata.write('\n') for item in constants_orig: if len(latex(item)) == 3: savedata.write(r'\def%s{\mathrm{%s}_\mathrm{%s}}' %(latex(item),latex(item)[1].lower(),latex(item)[2])) if len(latex(item)) == 4: savedata.write(r'\def%s{\mathrm{%s}_\mathrm{%s}}' %(latex(item),latex(item)[1].lower(),latex(item)[2:4])) savedata.write('\n') savedata.write('\n') for item in variables_dimless: if latex(item) == r"\tau": savedata.write(r'\def%s{\mathrm{%s}}' % (latex(item), latex(item)[1:4])) savedata.write('\n') else: savedata.write(r'\def%s{\mathrm{%s}}' %(latex(item),latex(item)[-1])) savedata.write('\n') for item in constants_dimless: if len(latex(item)) == 3: savedata.write(r'\def%s{\mathrm{%s}_\mathrm{%s}}' %(latex(item),latex(item)[1].lower(),latex(item)[2])) if len(latex(item)) == 4: savedata.write(r'\def%s{\mathrm{%s}_\mathrm{%s}}' %(latex(item),latex(item)[1].lower(),latex(item)[2:4])) savedata.write('\n') savedata.write('\chapter{Model 11}' + '\n') return def print_latex_file_ending(savedata,document_type="teildok"): if (document_type=="teildok"): savedata.write('\DocumentEnd' + '\n') savedata.close() else: savedata.write('\end{document}' + '\n') savedata.close() return #Think about sympify for a more flexible way to define symbols # First define the eqations as sting and sympify the string P,I,D,t = symbols("\P,\I,\D,\\t",nonnegative=True) Cp,Ci,Cd, Rp,Rpp,Rpi,Rpd, Ri,Rii,Rip,Rid, Rd,Rdd,Rdi = symbols("\Cp,\Ci,\Cd, \\Rp,\\Rpp,\\Rpi,\\Rpd \\Ri,\\Rii,\\Rip,\\Rid, \\Rd,\\Rdd,\\Rdi",nonnegative=True) variables_orig = [P,I,D,t] constants_orig = [Cp,Ci,Cd, Rp,Rpp,Rpi,Rpd, Ri,Rii,Rip,Rid, Rd,Rdd,Rdi] #print(variables_orig,constants_orig) # P -> X, I -> Y, D -> Z #dimensionless modell X,Y,Z,tau = symbols("\X,\Y,\Z,\\tau",nonnegative=True) #,positive=True cannot be specified, because there will be no X,Y,Z = 0 ... steady state found Rxy,Rxz, Ry,Ryx,Ryz, Rz,Rzy,Rzz = symbols("\\RXY,\\RXZ, \\RY,\\RYX,\\RYZ, \\RZ,\\RZY,\\RZZ",nonnegative=True) variables_dimless = [X,Y,Z,tau] constants_dimless = [Rxy,Rxz, Ry,Ryx,Ryz, Rz,Rzy,Rzz] #print(variables_dimless,constants_dimless) #start latex documentation document_name="test.tex" #document_type="teildok" document_type="standard_latex" savedata = open(document_name, "w") print_latex_header(savedata,variables_orig,constants_orig,variables_dimless,constants_dimless,document_type) #----------------------------------------------------------------------------------------------------------------------# # # # include an image of the system # # # #----------------------------------------------------------------------------------------------------------------------# #savedata.write('\section{System}'+'\n') #savedata.write('\\begin{figure}[!ht]'+'\n') #savedata.write('\\centering'+'\n') #savedata.write('\\includegraphics[width=0.4\\textwidth]{\string~/sysbio/arbeit/2016_05_02_modell1/bilder/modell11}'+ '\n') #savedata.write('\end{figure}'+'\n') #savedata.write('\\noindent\\rule{\\textwidth}{1pt}'+'\n') #----------------------------------------------------------------------------------------------------------------------# # # # DES of the system # # # #----------------------------------------------------------------------------------------------------------------------# dP = Rp * P * (1-P/Cp) dI = Ri * I * (1-I/Ci) + Rdi *D*I dD = Rd * D * (1-D/Cd) + Rpd *P*D + Rid *I*D savedata.write('\\begin{equation}'+'\n') savedata.write('\\begin{aligned}'+'\n') savedata.write('dP &= ' + latex(dP)+'\\\\ \n') savedata.write('dI &= ' + latex(dI)+'\\\\ \n') savedata.write('dD &= ' + latex(dD)+'\n') savedata.write('\end{aligned}'+'\n') savedata.write('\end{equation}'+'\n') savedata.write('\\noindent\\rule{\\textwidth}{1pt}'+'\n') #----------------------------------------------------------------------------------------------------------------------# # # # dimensionless DES # # # #----------------------------------------------------------------------------------------------------------------------# dimensionless_DES_list = [X*(1-X),Ry *(Y* (1 - Y) + Rzy * Z *Y ),Rz *(Z*(1-Z)+Rxz*X*Z+Ryz*Y*Z)] savedata.write('\\begin{equation}'+'\n') savedata.write('\\begin{aligned}'+'\n') for counter in range(len(dimensionless_DES_list)): savedata.write("d{} &= ".format(variables_dimless[counter]) + latex(dimensionless_DES_list[counter])+"\\\\ \n") savedata.write('\end{aligned}'+'\n') savedata.write('\end{equation}'+'\n') savedata.write('\\noindent\\rule{\\textwidth}{1pt}'+'\n') #----------------------------------------------------------------------------------------------------------------------# # # # parameter transformation # # # #----------------------------------------------------------------------------------------------------------------------# #X_subs = P /Cp #Y_subs = I /Ci #Z_subs = D /Cd #cp = rp/rpp #ci = ri/rii #cd = rd/rdd #Ry = Ri/Rp #Rzy = Rdi*Cd/Ri #Rz = Rd/Rp #Rxz = Rpd*Cp/Rd #Ryz = Rid*Ci/Rd #dP = dX.subs([( X, X_subs)]) #dI = dY.subs([( X, X_subs),( Y, Y_subs),( Z, Z_subs)]) #dI_output = (dI*Ci) #dI_output = dI_output.simplify() #print(dP) #print original model #----------------------------------------------------------------------------------------------------------------------# # # # compute the steady states and general Jacobian matrix # # # #----------------------------------------------------------------------------------------------------------------------# equilibria = solve( dimensionless_DES_list,(variables_dimless[0:len(dimensionless_DES_list)])) eqMat = Matrix(dimensionless_DES_list ) Mat = Matrix([variables_dimless[0:len(dimensionless_DES_list)]]) J = eqMat.jacobian(Mat) #print(J) savedata.write('\section{Steady States and Jacobian Matrix}'+'\n') savedata.write('\\begin{equation}'+'\n') savedata.write('\\begin{aligned}'+'\n') for counter in range(len(dimensionless_DES_list)): savedata.write('0 &= ' + latex(dimensionless_DES_list[counter].factor())+'\\\\ \n') savedata.write('\end{aligned}'+'\n') savedata.write('\end{equation}'+'\n') savedata.write('\\noindent\\rule{\\textwidth}{1pt}'+'\n') savedata.write('\\begin{equation}'+'\n') savedata.write('J = '+ latex(J)+'\n') savedata.write('\end{equation}'+'\n') savedata.write('\\noindent\\rule{\\textwidth}{1pt}'+'\n') for item in equilibria: J_sub = [(variables_dimless[counter],item[counter]) for counter in range(len(item)) ] eqmat = J.subs(J_sub) # substitute the equilibria points eqmat.simplify() #print(eqmat) J_sub = [] for counter in range(len(item)): if item[counter]!= 0 and item[counter]!= 1: J_sub.append((item[counter],variables_dimless[counter])) for i in range(len(item)): print(eqmat[i,i])# = eqmat[i][i].subs((item[counter],-1* variables_dimless[counter])) print(solve_univariate_inequality(eqmat[i,i] >= 0, (Rxy,Rxz, Ry,Ryx,Ryz, Rz,Rzy,Rzz))) eqmat = eqmat.subs(J_sub) #print(J_sub) #print(eqmat) savestring = '\\textbf{steady state} $(' for counter in range(len(item)): savestring += latex(item[counter]) + ',' savestring += ')$'+ '\n' savedata.write(savestring) savedata.write('\\begin{equation}' + '\n') savedata.write('J = ' + latex(eqmat) + '\n') savedata.write('\end{equation}' + '\n') savedata.write('\\begin{equation}' + '\n') savedata.write('\\begin{aligned}' + '\n') for counter in range(len(item)): savedata.write('\lambda_{} &= '.format(counter+1) + latex(eqmat.eigenvals().keys()[counter]) + '\\\\ \n') savedata.write('\end{aligned}' + '\n') savedata.write('\end{equation}' + '\n') savedata.write('\\noindent\\rule{\\textwidth}{1pt}'+'\n') print('The eigenvalues for the fixed point (%s, %s, %s) are %s, %s and %s:' %(item[0], item[1], item[2], eqmat.eigenvals().keys()[0], eqmat.eigenvals().keys()[1], eqmat.eigenvals().keys()[2])) print('-------------------------------------------') print_latex_file_ending(savedata,document_type)
#!/bin/bash ./test.py pdflatex test.tex evince test.pdf &
