Ask Your Question

Substituting a differential equation into an expression

asked 2018-11-01 07:57:04 +0200

Hans gravatar image

updated 2018-11-01 14:44:37 +0200

Say I have a (heat) partial differential equation of $f(t,x)$ $$\frac{\partial f}{\partial t}=\frac12\frac{\partial^2 f}{\partial x^2}.$$ I would like to express an algebraic function $g$ of the variable $\Big(\frac{\partial f}{\partial t}, \frac{\partial^2 f}{\partial t\partial x}\Big)$ all in the partial derivatives of $x$. That is $$\frac{\partial^2 f}{\partial t\partial x}=\frac12\frac{\partial^3 f}{\partial x^3}.$$ The result should be $$g\Big(\frac{\partial f}{\partial t}, \frac{\partial^2 f}{\partial t\partial x}\Big)=g\Big(\frac12\frac{\partial^2 f}{\partial x^2}, \frac12\frac{\partial^3 f}{\partial x^3}\Big).$$

How can I accomplish this with SageMath?

edit retag flag offensive close merge delete

1 Answer

Sort by ยป oldest newest most voted

answered 2018-11-01 17:28:15 +0200

rburing gravatar image

updated 2020-02-14 23:18:57 +0200

This is tricky, but it can be done. The idea is to replace (in the expression tree) derivatives of $f$ (with at least one derivative with respect to $t$) by the right-hand side of the equation, substituted into itself as many times as there are derivatives with respect to $t$ minus one, and finally replacing $f$ by the derivative of $f$ with respect to the remaining variables (if any).

My attempt to implement this (based on SubstituteFunction) is below, with an example (your example) in the docstring.

from sage.symbolic.expression_conversions import ExpressionTreeWalker
from sage.symbolic.operators import FDerivativeOperator

class SubstituteEvolutionaryPDE(ExpressionTreeWalker):
    def __init__(self, pde, t):
        A class that walks the tree and replaces derivatives of `f` with
        respect to `t` by the right-hand side of an evolutionary PDE.
            sage: var('x,t'); f = function('f'); g = function('g')
            sage: pde = diff(f(t,x),t) == 1/2*diff(f(t,x),x,x)
            sage: s = SubstituteEvolutionaryPDE(pde, t)
            sage: h = g(diff(f(t,x),t), diff(f(t,x),t,x))
            sage: s(h)
            g(1/2*diff(f(t, x), x, x), 1/2*diff(f(t, x), x, x, x))


        ``pde`` is of the form ``diff(f(t,...),t) == ...`` and the first
        argument of ``f`` is always ``t``.
        f_expr = integrate(pde.lhs(), t)
        self.f = f_expr.operator()
        self.args = f_expr.operands()
        self.rhs = pde.rhs()

    def derivative(self, ex, operator):
        if operator.function() == self.f and 0 in operator.parameter_set(): # f is differentiated with respect to t
            t_derivatives = [p for p in operator.parameter_set() if p == 0] # (assumes t is the first argument)
            result = self.rhs
            for _ in range(len(t_derivatives)-1):
                result = result.substitute_function(self.f, self.rhs.function(*self.args))
            if len(t_derivatives) < len(operator.parameter_set()): # derivatives w.r.t variables other than t
                other_derivatives = [p for p in operator.parameter_set() if not p == 0] # (assumes t is the first argument)
                new_operator = FDerivativeOperator(operator.function(), other_derivatives)
                result = result.substitute_function(self.f, new_operator(*[self(_) for _ in ex.operands()]).function(*self.args))
            return result
            return operator(*[self(_) for _ in ex.operands()])

Edit: updated to use substitute_function correctly (the second argument should also be a function).

edit flag offensive delete link more


Impressive! I am digesting your answer. What does the * in 'new_operator(*[self(_) for _ in ex.operands()])' do?

Hans gravatar imageHans ( 2018-11-01 20:57:53 +0200 )edit

This "unpacks" the list so e.g. f(*[t,x]) == f(t,x); it's a way to insert positional arguments from a (possibly dynamically generated, e.g. of variable length) list.

rburing gravatar imagerburing ( 2018-11-01 21:21:12 +0200 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower


Asked: 2018-11-01 07:57:04 +0200

Seen: 543 times

Last updated: Feb 14 '20