Hello,
I am currently converting, for teaching purposes, a number of code examples
from Mathematica to sympy in order to have a purely Python-based teaching
environment (the numerical part has long been done in Numpy/Scipy).
I have now run into a situation where I can get some code to work in a
mathematically correct way, but produce output that makes me think I am not
doing things quite right.
The task is to use sympy to verify the order of a finite difference time
integrator for an ODE (for the sake of simplicity, autonomous), which
involves Taylor expanding the error with respect to the stepsize parameter,
the substituting the differential equation and its derivatives.
The attached code does just that, here with the implicit midpoint rule as a
simple example. The issue I have is that sympy does not seem to have the
notion of an abstract derivative, so the series expansion produces terms
where Subs is used to represent derivatives with an argument that is not
simply a symbol. In the end, however, I have to force evaluation with
doit((), so that terms cancel properly, which has the side effect that all
Subs are resolved and I get nonsense terms like
Derivative(f(y(0)), y(0))
Of course, I could just ignore this as I have already achieved my goal, but
I am trying to teach "good" sympy and this result strikes me as not being
formulated the way it should.
Any comments highly appreciated!
Marcel
--
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 view this discussion on the web visit
https://groups.google.com/d/msgid/sympy/fd11673e-c17a-4131-aee4-38ab2ea00661n%40googlegroups.com.
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Oct 5 14:35:11 2021
@author: oliver
"""
from sympy import *
init_printing()
h,t = symbols("h,t")
f = Function("f")
y = Function("y")
a = Wild('a')
# Set up the local error expression for the midpoint scheme
# for an autonomous ODE
error = y(0) + h*f((y(0)+y(h))/2) - y(h)
# Now Taylor expand with respect to h:
p = 3 # Compute explicitly up to order p
error = series(error, h, n=p+1)
# Now recursively substitute derivatives of the differential equation,
# starting with the highest possible derivative and ending with the
# differential equation itself:
for i in range(p,1,-1):
error = error.replace(Derivative(y(a),(a,i)), Subs(diff(f(y(t)),(t,i-1)), t, a))
error = error.replace(Derivative(y(a),a),f(y(a)))
print(error.doit())
