# ## MIT License

# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
#

# In[1]:
import sympy as s
s.init_printing()
# Establish some variables that we'll need:
# In[6]:
u = s.Function("u")
a, x, t, h_x, h_t = s.symbols("a, x, t, h_x, h_t")
xi, xi_1, xi_2, tau = s.symbols("xi, xi1, xi2, tau")
# `taylor` is a utility function that spits out a taylor expansion for $f(x+h)$, optionally including a remainder term, with all variables under our control.
# In[7]:
def taylor(f, x, h, n, remainder_variable=None):
result = sum(f.diff(x, i)*h**i/s.factorial(i) for i in range(n))
if remainder_variable:
result += f.diff(x, n).subs(x, remainder_variable)*h**n/s.factorial(n)
return result
# - Try it out by expanding $u(x+h_x,t)$
# - Vary the order
# - Expand $u(x,t+h_t)$ instead
# In[9]:
taylor(u(x,t), x, h_x, 3, xi)
# Assign the PDE we're solving to `pde`:
# In[10]:
pde = u(x, t).diff(t) + a * u(x, t).diff(x)
pde
# Write out the scheme we're analyzing, in this case ETCS:
# In[11]:
etcs = (
(u(x, t+h_t) - u(x, t))/h_t
+
a*(u(x+h_x, t) - u(x-h_x, t))/(2*h_x))
etcs
# Follow this general pattern:
# ```
# etcs
# .subs(u(x, t+h_t), taylor(u(x,t), t, h_t, 2, tau))
# ```
# to arrive at the truncation error.
#
# ⚠️ Make sure to keep the two $x$ remainder terms separate.
# In[13]:
etcs_taylor = (
etcs
.subs(u(x, t+h_t), taylor(u(x,t), t, h_t, 2, tau))
.subs(u(x+h_x, t), taylor(u(x,t), x, h_x, 3, xi_1))
.subs(u(x-h_x, t), taylor(u(x,t), x, -h_x, 3, xi_2))
)
s.simplify(etcs_taylor - pde)
# In[ ]: