Explicitly solved backward Euler step

backward_euler.ipynb

{"nbformat_minor": 0, "cells": [{"execution_count": 1, "cell_type": "code", "source": "from sympy import *", "outputs": [], "metadata": {"collapsed": true, "trusted": true}}, {"execution_count": 5, "cell_type": "code", "source": "'''\ndC1/dt = Gr0*C22 + Gd1*O1 - Ga1*C1 : 1\ndO1/dt = Ga1*C1 - (Gd1+Gf)*O1 + Gb*O2 : 1\ndO2/dt = Ga2*C22 + Gf*O1 - (Gd2+Gb)*O2 : 1\nC22 = 1 - O1 - O2 - C1 : 1\n'''\nC1, O1, O2 = symbols('C1 O1 O2', real=True)\nC1_prev, O1_prev, O2_prev = symbols('C1_prev O1_prev O2_prev', real=True)\nGr0, Gd1, Gd2, Ga1, Ga2, Gb, Gf = symbols('Gr0 Gd1 Gd2 Ga1 Ga2 Gb Gf', real=True)", "outputs": [], "metadata": {"collapsed": true, "trusted": true}}, {"execution_count": 6, "cell_type": "code", "source": "# Differential equations with C22 substituted\ndiff_eq1 = Gr0*(1 - O1 - O2 - C1) + Gd1*O1 - Ga1*C1\ndiff_eq2 = Ga1*C1 - (Gd1+Gf)*O1 + Gb*O2\ndiff_eq3 = Ga2*(1 - O1 - O2 - C1) + Gf*O1 - (Gd2+Gb)*O2", "outputs": [], "metadata": {"collapsed": true, "trusted": true}}, {"execution_count": 7, "cell_type": "code", "source": "dt = Symbol('dt', real=True, positive=True)\neq1 = Eq(C1, C1_prev + dt * diff_eq1)\neq2 = Eq(O1, O1_prev + dt * diff_eq2)\neq3 = Eq(O2, O2_prev + dt * diff_eq3)\n", "outputs": [], "metadata": {"collapsed": true, "trusted": true}}, {"execution_count": 9, "cell_type": "code", "source": "solution = solve([eq1, eq2, eq3], C1, O1, O2)", "outputs": [], "metadata": {"collapsed": true, "trusted": true}}, {"execution_count": 10, "cell_type": "code", "source": "solution", "outputs": [{"execution_count": 10, "output_type": "execute_result", "data": {"text/plain": "{C1: -(O1_prev*dt*(Gr0*dt*(Ga2 - Gf) + (Gd1 - Gr0)*(Ga2*dt + Gb*dt + Gd2*dt + 1)) + dt*(Ga2*dt + O2_prev)*(Gb*dt*(Gd1 - Gr0) - Gr0*(Gd1*dt + Gf*dt + 1)) + (C1_prev + Gr0*dt)*(Gb*dt**2*(Ga2 - Gf) + (Gd1*dt + Gf*dt + 1)*(Ga2*dt + Gb*dt + Gd2*dt + 1)))/(Ga1*Gr0*dt**3*(Ga2 - Gf) + Ga1*dt**2*(Gd1 - Gr0)*(Ga2*dt + Gb*dt + Gd2*dt + 1) - Ga2*Gb*dt**3*(Gd1 - Gr0) + Ga2*Gr0*dt**2*(Gd1*dt + Gf*dt + 1) - Gb*dt**2*(Ga2 - Gf)*(Ga1*dt + Gr0*dt + 1) - (Ga1*dt + Gr0*dt + 1)*(Gd1*dt + Gf*dt + 1)*(Ga2*dt + Gb*dt + Gd2*dt + 1)),\n O2: (O1_prev*dt*(Ga2*dt*(Gd1 - Gr0) + (Ga2 - Gf)*(Ga1*dt + Gr0*dt + 1)) + dt*(C1_prev + Gr0*dt)*(Ga1*dt*(Ga2 - Gf) + Ga2*(Gd1*dt + Gf*dt + 1)) + (Ga2*dt + O2_prev)*(Ga1*dt**2*(Gd1 - Gr0) - (Ga1*dt + Gr0*dt + 1)*(Gd1*dt + Gf*dt + 1)))/(Ga1*Gr0*dt**3*(Ga2 - Gf) + Ga1*dt**2*(Gd1 - Gr0)*(Ga2*dt + Gb*dt + Gd2*dt + 1) - Ga2*Gb*dt**3*(Gd1 - Gr0) + Ga2*Gr0*dt**2*(Gd1*dt + Gf*dt + 1) - Gb*dt**2*(Ga2 - Gf)*(Ga1*dt + Gr0*dt + 1) - (Ga1*dt + Gr0*dt + 1)*(Gd1*dt + Gf*dt + 1)*(Ga2*dt + Gb*dt + Gd2*dt + 1)),\n O1: (O1_prev*(Ga2*Gr0*dt**2 - (Ga1*dt + Gr0*dt + 1)*(Ga2*dt + Gb*dt + Gd2*dt + 1)) - dt*(C1_prev + Gr0*dt)*(Ga1*(Ga2*dt + Gb*dt + Gd2*dt + 1) - Ga2*Gb*dt) + dt*(Ga2*dt + O2_prev)*(Ga1*Gr0*dt - Gb*(Ga1*dt + Gr0*dt + 1)))/(Ga1*Gr0*dt**3*(Ga2 - Gf) + Ga1*dt**2*(Gd1 - Gr0)*(Ga2*dt + Gb*dt + Gd2*dt + 1) - Ga2*Gb*dt**3*(Gd1 - Gr0) + Ga2*Gr0*dt**2*(Gd1*dt + Gf*dt + 1) - Gb*dt**2*(Ga2 - Gf)*(Ga1*dt + Gr0*dt + 1) - (Ga1*dt + Gr0*dt + 1)*(Gd1*dt + Gf*dt + 1)*(Ga2*dt + Gb*dt + Gd2*dt + 1))}"}, "metadata": {}}], "metadata": {"collapsed": false, "trusted": true}}, {"execution_count": 15, "cell_type": "code", "source": "from brian2.parsing.sympytools import sympy_to_str\nfrom brian2.utils.stringtools import word_substitute", "outputs": [], "metadata": {"collapsed": false, "trusted": true}}, {"execution_count": 23, "cell_type": "code", "source": "# Generate abstract code for Brian 2 \nreplacements = {'C1_prev': 'C1', 'O1_prev': 'O1', 'O2_prev': 'O2'}\nprint '_C1 = ' + word_substitute(sympy_to_str(solution[C1]), replacements)\nprint '_O1 = ' + word_substitute(sympy_to_str(solution[O1]), replacements)\nprint '_O2 = ' + word_substitute(sympy_to_str(solution[O2]), replacements)\nprint 'C1 = _C1'\nprint 'O1 = _O1'\nprint 'O2 = _O2'", "outputs": [{"output_type": "stream", "name": "stdout", "text": "_C1 = -(O1*dt*(Gr0*dt*(Ga2 - Gf) + (Gd1 - Gr0)*(Ga2*dt + Gb*dt + Gd2*dt + 1)) + dt*(Ga2*dt + O2)*(Gb*dt*(Gd1 - Gr0) - Gr0*(Gd1*dt + Gf*dt + 1)) + (C1 + Gr0*dt)*(Gb*dt**2*(Ga2 - Gf) + (Gd1*dt + Gf*dt + 1)*(Ga2*dt + Gb*dt + Gd2*dt + 1)))/(Ga1*Gr0*dt**3*(Ga2 - Gf) + Ga1*dt**2*(Gd1 - Gr0)*(Ga2*dt + Gb*dt + Gd2*dt + 1) - Ga2*Gb*dt**3*(Gd1 - Gr0) + Ga2*Gr0*dt**2*(Gd1*dt + Gf*dt + 1) - Gb*dt**2*(Ga2 - Gf)*(Ga1*dt + Gr0*dt + 1) - (Ga1*dt + Gr0*dt + 1)*(Gd1*dt + Gf*dt + 1)*(Ga2*dt + Gb*dt + Gd2*dt + 1))\n_O1 = (O1*(Ga2*Gr0*dt**2 - (Ga1*dt + Gr0*dt + 1)*(Ga2*dt + Gb*dt + Gd2*dt + 1)) - dt*(C1 + Gr0*dt)*(Ga1*(Ga2*dt + Gb*dt + Gd2*dt + 1) - Ga2*Gb*dt) + dt*(Ga2*dt + O2)*(Ga1*Gr0*dt - Gb*(Ga1*dt + Gr0*dt + 1)))/(Ga1*Gr0*dt**3*(Ga2 - Gf) + Ga1*dt**2*(Gd1 - Gr0)*(Ga2*dt + Gb*dt + Gd2*dt + 1) - Ga2*Gb*dt**3*(Gd1 - Gr0) + Ga2*Gr0*dt**2*(Gd1*dt + Gf*dt + 1) - Gb*dt**2*(Ga2 - Gf)*(Ga1*dt + Gr0*dt + 1) - (Ga1*dt + Gr0*dt + 1)*(Gd1*dt + Gf*dt + 1)*(Ga2*dt + Gb*dt + Gd2*dt + 1))\n_O2 = (O1*dt*(Ga2*dt*(Gd1 - Gr0) + (Ga2 - Gf)*(Ga1*dt + Gr0*dt + 1)) + dt*(C1 + Gr0*dt)*(Ga1*dt*(Ga2 - Gf) + Ga2*(Gd1*dt + Gf*dt + 1)) + (Ga2*dt + O2)*(Ga1*dt**2*(Gd1 - Gr0) - (Ga1*dt + Gr0*dt + 1)*(Gd1*dt + Gf*dt + 1)))/(Ga1*Gr0*dt**3*(Ga2 - Gf) + Ga1*dt**2*(Gd1 - Gr0)*(Ga2*dt + Gb*dt + Gd2*dt + 1) - Ga2*Gb*dt**3*(Gd1 - Gr0) + Ga2*Gr0*dt**2*(Gd1*dt + Gf*dt + 1) - Gb*dt**2*(Ga2 - Gf)*(Ga1*dt + Gr0*dt + 1) - (Ga1*dt + Gr0*dt + 1)*(Gd1*dt + Gf*dt + 1)*(Ga2*dt + Gb*dt + Gd2*dt + 1))\nC1 = _C1\nO1 = _O1\nO2 = _O2\n"}], "metadata": {"collapsed": false, "trusted": true}}, {"execution_count": null, "cell_type": "code", "source": "", "outputs": [], "metadata": {"collapsed": true, "trusted": true}}], "nbformat": 4, "metadata": {"kernelspec": {"display_name": "Python 2", "name": "python2", "language": "python"}, "language_info": {"mimetype": "text/x-python", "nbconvert_exporter": "python", "version": "2.7.10", "name": "python", "file_extension": ".py", "pygments_lexer": "ipython2", "codemirror_mode": {"version": 2, "name": "ipython"}}}}