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 sympy+unsubscr...@googlegroups.com. 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())