Starting from input parameters `sigma0`, `sigmadot0`, `sigmaddot0`, `psi2`, `psi1`, and `psi0`, representing the corresponding values of the Bondi data at time $u=0$, we want to integrate the initial-value problem to find all the Bondi data at all times.  We start with
\begin{equation}
  \sigma = \texttt{sigma0} + u\,\texttt{sigmadot0} + u^2\,\texttt{sigmaddot0}\,/\,2,
\end{equation}
impose an initial condition on $\psi_2$ given by
\begin{equation}
  \Im[ψ₂] = -\Im[ð²σ̄ + σ ∂σ̄/∂u]
\end{equation}
so that the mass parameter is real, and then solve for the remaining parameters according to
\begin{align}
    ψ̇₀ &= -ðψ₁ + 3 σ ψ₂ \\
    ψ̇₁ &= -ðψ₂ + 2 σ ψ₃ \\
    ψ̇₂ &= -ðψ₃ + 1 σ ψ₄ \\
    ψ₃ &= ð ∂σ̄/∂u \\
    ψ₄ &= - ∂²σ̄/∂u²
\end{align}

Sympy can do all the heavy lifting for us, including writing the actual code.  There's a small difficulty with getting $\eth$ to act like a linear operator, which is why we have the Ð function below.


```python
import sympy
from sympy import init_printing
init_printing()

u = sympy.symbols('u', real=True)
ð = sympy.symbols('ð', cls=sympy.Function)

dot = lambda f: sympy.Derivative(f, u, evaluate=True)
ddot = lambda f: sympy.Derivative(f, u, u, evaluate=True)
integral = lambda f: sympy.Integral(f, u).doit()
```


```python
def Ð(expr, expand=True):
    """Apply ð and expand by linearity over numbers and u**n
    
    This is inspired by https://stackoverflow.com/a/56924562/1194883,
    which makes clever use of sympy's `expr.replace(filter, func)`
    functionality to expand sums inside the argument of ð, and then
    I make dumb use of the basic replace functionality (with
    wildcards) to expand scalar multiplication, assuming the scalars
    are just literal numbers and powers of `u`.
    
    """
    import sympy
    c = sympy.Wild('c', properties=(lambda x: x.is_Number,))
    p = sympy.Wild('p')
    w = sympy.Wild('w')
    expr = ð(expr)
    return expr.expand().replace(
        lambda x: x.func == ð,
        lambda x: sympy.Add(*[x.func(a) for a in sympy.Add.make_args(x.args[0])])
    ).replace(ð(c * u**p * w), c * u**p * ð(w))

```


```python
def split_polynomial(expr, string):
    import sympy
    from sympy.polys.polyfuncs import horner
    from sympy.printing.pycode import pycode
    from sympy.codegen.ast import Assignment
    expr = sympy.Poly(expr, u)
    values = [sympy.factor(c) for c in reversed(expr.all_coeffs())]
    n_coeffs = len(values)
    symbols = sympy.symbols(string+f"_:{n_coeffs}")
    symbol = sympy.symbols(string)
    value = horner(sympy.Poly(reversed(symbols), u))
    print('\n'.join(pycode(Assignment(s, v), allow_unknown_functions=True) for s, v in zip(symbols, values)))
    print(pycode(Assignment(symbol, value)))
    return value
```


```python
print('ð = lambda x: x.eth')
print('conjugate = lambda x: x.bar')
σ = split_polynomial(sympy.symbols('sigma0') + u * sympy.symbols('sigmadot0') + u**2 * sympy.symbols('sigmaddot0') / 2, 'σ')
print('# ψ₄ = -∂²σ̄/∂u²')
ψ4 = split_polynomial(-ddot(σ.conjugate()), 'ψ4')
print('# ψ₃ = ð ∂σ̄/∂u')
ψ3 = split_polynomial(Ð(dot(σ.conjugate())), 'ψ3')
print('# ψ₂ = ∫ (-ðψ₃ + σψ₄) du')
ψ̇2 = sympy.collect(sympy.expand(-Ð(ψ3) + σ*ψ4), u)
ψ2 = split_polynomial(integral(ψ̇2) + sympy.symbols('psi2'), 'ψ2')
print('# ψ₁ = ∫ (-ðψ₂ + 2σψ₃) du')
ψ̇1 = sympy.collect(sympy.expand(-Ð(ψ2) + 2*σ*ψ3), u)
ψ1 = split_polynomial(integral(ψ̇1) + sympy.symbols('psi1'), 'ψ1')
print('# ψ₀ = ∫ (-ðψ₁ + 3σψ₂) du')
ψ̇0 = sympy.collect(sympy.expand(-Ð(ψ1) + 3*σ*ψ2), u)
ψ0 = split_polynomial(integral(ψ̇0) + sympy.symbols('psi0'), 'ψ0')
```

    ð = lambda x: x.eth
    conjugate = lambda x: x.bar
    σ_0 = sigma0
    σ_1 = sigmadot0
    σ_2 = (1/2)*sigmaddot0
    σ = u*(u*σ_2 + σ_1) + σ_0
    # ψ₄ = -∂²σ̄/∂u²
    ψ4_0 = -2*conjugate(σ_2)
    ψ4 = ψ4_0
    # ψ₃ = ð ∂σ̄/∂u
    ψ3_0 = ð(conjugate(σ_1))
    ψ3_1 = 2*ð(conjugate(σ_2))
    ψ3 = u*ψ3_1 + ψ3_0
    # ψ₂ = ∫ (-ðψ₃ + σψ₄) du
    ψ2_0 = psi2
    ψ2_1 = σ_0*ψ4_0 - ð(ψ3_0)
    ψ2_2 = -1/2*(-σ_1*ψ4_0 + ð(ψ3_1))
    ψ2_3 = (1/3)*σ_2*ψ4_0
    ψ2 = u*(u*(u*ψ2_3 + ψ2_2) + ψ2_1) + ψ2_0
    # ψ₁ = ∫ (-ðψ₂ + 2σψ₃) du
    ψ1_0 = psi1
    ψ1_1 = 2*σ_0*ψ3_0 - ð(ψ2_0)
    ψ1_2 = -1/2*(-2*σ_0*ψ3_1 - 2*σ_1*ψ3_0 + ð(ψ2_1))
    ψ1_3 = -1/3*(-2*σ_1*ψ3_1 - 2*σ_2*ψ3_0 + ð(ψ2_2))
    ψ1_4 = -1/4*(-2*σ_2*ψ3_1 + ð(ψ2_3))
    ψ1 = u*(u*(u*(u*ψ1_4 + ψ1_3) + ψ1_2) + ψ1_1) + ψ1_0
    # ψ₀ = ∫ (-ðψ₁ + 3σψ₂) du
    ψ0_0 = psi0
    ψ0_1 = 3*σ_0*ψ2_0 - ð(ψ1_0)
    ψ0_2 = -1/2*(-3*σ_0*ψ2_1 - 3*σ_1*ψ2_0 + ð(ψ1_1))
    ψ0_3 = -1/3*(-3*σ_0*ψ2_2 - 3*σ_1*ψ2_1 - 3*σ_2*ψ2_0 + ð(ψ1_2))
    ψ0_4 = -1/4*(-3*σ_0*ψ2_3 - 3*σ_1*ψ2_2 - 3*σ_2*ψ2_1 + ð(ψ1_3))
    ψ0_5 = -1/5*(-3*σ_1*ψ2_3 - 3*σ_2*ψ2_2 + ð(ψ1_4))
    ψ0_6 = (1/2)*σ_2*ψ2_3
    ψ0 = u*(u*(u*(u*(u*(u*ψ0_6 + ψ0_5) + ψ0_4) + ψ0_3) + ψ0_2) + ψ0_1) + ψ0_0



```python

```

d6851f65cfcf97394ef5b754e1757ed220c2e642
ipynb
Jupyter Notebook
notes/bondi_initial_value_solution.ipynb
akhairna/scri
3b7f307d19ef303914cef2fa088ee750ef8533c2
MIT
2018-03-11T03:09:19.000Z
2022-03-08T20:09:25.000Z
notes/bondi_initial_value_solution.ipynb
akhairna/scri
3b7f307d19ef303914cef2fa088ee750ef8533c2
MIT
2016-08-11T17:31:12.000Z
2021-11-17T17:18:23.000Z
notes/bondi_initial_value_solution.ipynb
akhairna/scri
3b7f307d19ef303914cef2fa088ee750ef8533c2
MIT
2015-06-22T22:16:30.000Z
2022-02-22T23:38:26.000Z
Qwen/Qwen-72B
1. YES
2. YES