#!/usr/bin/env python
# coding: utf-8

# # BVP Green's Functions
# 
# Copyright (C) 2026 Andreas Kloeckner
# 
# <details>
# <summary>MIT License</summary>
# 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.
# </details>

# In[2]:


import numpy as np
import numpy.linalg as la
import sympy as sp
import matplotlib.pyplot as plt


# In[3]:


x = sp.Symbol("x")
z = sp.Symbol("z")


# Let's consider the ODE system
# $$
# \boldsymbol y'(x) = \begin{bmatrix}
# 0 & 1\\
# -1 & 0
# \end{bmatrix}
# \boldsymbol y(x) + \begin{bmatrix}
# 3\\
# 1
# \end{bmatrix}
# $$
# i.e. the harmonic oscillator with boundary conditions
# $$
# \begin{align*}
# y_1(0) &= 2\\
# y_1'(\pi/2) - 3 y_1(\pi/2) &= -4\\
# \end{align*}
# $$

# Set up system data and the fundamental solution matrix $\boldsymbol Y$:

# In[4]:


A = sp.Matrix([[0, 1], [-1, 0]])
rhs = sp.Matrix([[3, 1]]).T
a = sp.sympify(0)
b = sp.pi/2
B_a = sp.Matrix([[1, 0], [0, 0]])
B_b = sp.Matrix([[0, 0], [-3, 1]])
c = sp.Matrix([[2, -4]]).T

xmesh = np.linspace(float(a), float(b), 1000)

solutions = [sp.sin(x), sp.cos(x)]
Y = sp.Matrix([
   [sol.diff(x, num_derivs) for sol in solutions]
   for num_derivs in [0, 1]
])
Y


# Check that the columns satisfies the homogeneous ODE:

# In[5]:


Y.diff(x) - A @ Y


# ## Solving the volume-homogeneous/"boundary-only" BVP 
# 
# Plug in the BCs. 

# In[6]:


B_a @ Y.subs(x, a) 


# In[7]:


B_b @ Y.subs(x, b)


# What do the matrix entries of $\boldsymbol Q$ mean?

# In[8]:


Q = B_a @ Y.subs(x, a) + B_b @ Y.subs(x, b)
Q


# What does $\boldsymbol Q^{-1} \boldsymbol c$ mean?

# In[9]:


Q.inv() @ c


# In[10]:


Phi = Y @ Q.inv()
Phi


# Observe:

# In[112]:


B_a @ Phi.subs(x, a) + B_b @ Phi.subs(x, b)


# Specifically:

# In[115]:


B_a @ Phi.subs(x, a)


# In[114]:


B_b @ Phi.subs(x, b)


# Check $\boldsymbol y_B = \boldsymbol \Phi \boldsymbol c$ satisfies the BCs:

# In[11]:


yB = Phi @ c
yB


# In[12]:


B_a @ yB.subs(x, a)


# In[13]:


B_b @ yB.subs(x, b)


# Plot it:

# In[14]:


plt.plot(xmesh, sp.lambdify(x, yB[0], "numpy")(xmesh))


# ## Solving the boundary-homogeneous/"volume-only" BVP 

# Set up the two pieces of the Green's function:

# In[104]:


Gleft = Phi @ B_a @ Phi.subs(x, a) @ Phi.subs(x, z).inv()
Gright = - Phi @ B_b @ Phi.subs(x, b) @ Phi.subs(x, z).inv()
Gleft


# Define an overall function and plot:

# In[105]:


G0 = sp.Piecewise(
    (Gleft[0, 0], z < x),
    (Gright[0, 0], z > x),
)
plt.plot(xmesh, sp.lambdify(x, G0.subs(z, 0.9), "numpy")(xmesh))
G0


# Verify that G satisfies the volume-homogeneous ODE:

# In[106]:


sp.simplify(Gleft.diff(x) - A @ Gleft)


# In[107]:


sp.simplify(Gright.diff(x) - A @ Gright)


# Verify that G satisfies the boundary conditions.
# 
# - Note that $x= a<z$, so that we use the $x < z$ `Gright` branch.
# - 

# In[108]:


B_a @ Gright.subs(x, a)


# In[109]:


B_b @ Gleft.subs(x, b)


# Verify the jump condition:

# In[101]:


sp.simplify(Gleft.subs(x, z) - Gright.subs(x, z))


# Put together the volume solution:

# In[86]:


yV =  (
    sp.integrate(Gleft @ rhs.subs(x, z), (z, a, x))
    + sp.integrate(Gright @ rhs.subs(x, z), (z, x, b))
)
yV


# Verify the ODE:

# In[87]:


sp.simplify(yV.diff(x) - A @ yV - rhs)


# OK then...
# Verify the ODE numerically:

# In[88]:


sp.simplify(yV.diff(x) - A @ yV - rhs).subs(x, 0.2).evalf()


# Plot it:

# In[89]:


plt.plot(xmesh, sp.lambdify(x, yV[0], "numpy")(xmesh))


#