This comes from a Hamilton-Jacobi-Bellman equation and I cannot isolate the second derivative, so I guess that makes it a differential algebraic equation. I have been trying to use assimulo but I can't manage to solve for Psi. The idea is to use the shooting method by choosing the boundary $\bar{\omega}$ correctly and trying to match Psi(0)=L. Does anyone know how to do this? What ressources I could use to properly code this? import numpy as np import matplotlib.pyplot as plt from scipy.interpolate import CubicSpline from scipy.integrate import solve_ivp from scipy.optimize import fsolve from assimulo.solvers import IDA from assimulo.problem import Implicit_Problem from assimulo.exception import TerminateSimulation # === Parameters === bar_omega = 0.6 L = 0.3 r = 0.045 mu = 0.08 gamma = 0.048 sigma_X = 1 sigma_Y = 1 rho = 0.5 alpha = 0.25 phi = 0.2 lambda_theta = 1 lambda_nu = 0.25 xi = 0.1 thetamin = 0 thetamax = xi / alpha numin = 0 numax = r / alpha # === Helper functions === def nu_star(omega, psi, dpsi, ddpsi): t1 = mu * (1 + psi - omega * dpsi) t2 = rho * sigma_X * sigma_Y * lambda_nu * lambda_theta * ddpsi t3 = alpha * phi + rho * sigma_X * sigma_Y * lambda_theta * omega * ddpsi t4 = alpha * lambda_theta - (lambda_theta**2) * (sigma_X**2) * ddpsi t5 = omega * lambda_nu * (sigma_Y**2) * ddpsi numerator = t1 - t2 * (t3 / t4) - t5 a1 = mu * lambda_nu - (lambda_nu**2) * (sigma_Y**2) * ddpsi a2 = (t2**2) / t4 denominator = a1 - a2 return numerator / denominator def theta_star(omega, psi, dpsi, ddpsi): nu = nu_star(omega, psi, dpsi, ddpsi) numerator = alpha * phi - rho * sigma_X * sigma_Y * lambda_theta * (lambda_nu * nu - omega) * ddpsi denominator = alpha * lambda_theta - (lambda_theta**2) * (sigma_X**2) * ddpsi return numerator / denominator def HJB_RES(omega, psi_val, dpsi_val, ddpsi): psi = lambda o: psi_val dpsi = lambda o: dpsi_val nu = nu_star(omega, psi_val, dpsi_val, ddpsi) theta = theta_star(omega, psi_val, dpsi_val, ddpsi) term1 = (r - mu * nu) * psi(omega) term2 = omega * (gamma - mu * nu) * dpsi(omega) vol = ((lambda_theta * theta * sigma_X)**2 + sigma_Y**2 * (lambda_nu * nu - omega)**2 - 2 * rho * sigma_X * sigma_Y * lambda_theta * theta * (lambda_nu * nu - omega)) term3 = 0.5 * vol * ddpsi term4 = (-mu * nu + phi * (xi - alpha * theta) + 0.5 * (lambda_nu * mu * nu**2 + lambda_theta * alpha * theta**2)) return term1 - term2 - term3 + term4 def bar_Psi(bar_omega): b1 = (r * lambda_nu - mu * (mu + bar_omega)) / mu b2 = np.sqrt(((lambda_nu / mu) * r + 1 - mu)**2 + 2 * (lambda_nu / mu) * (bar_omega * (gamma - r) - r + phi * xi - 0.5 * alpha * (phi**2 / lambda_theta))) return b1 + b2 # === Define DAE system === def equations(t, y, yprime): print(f"Residual eval at t={t:.6f}") psi_val = y[0] dpsi_val = y[1] dpsi_t = yprime[0] ddpsi = yprime[1] omega = t res0 = dpsi_val - dpsi_t res1 = HJB_RES(omega, psi_val, dpsi_val, ddpsi) return np.array([res0, res1]) # === Initial/Final conditions === y_f = [bar_Psi(bar_omega), -1] yprime_f = [-1, 0] t0 = bar_omega # Set up the problem and solver model = Implicit_Problem(equations, y_f, yprime_f, t0) sim = IDA(model) sim.backward = True sim.algvar = [True, True] sim.rtol = 1e-6 sim.atol = [1e-6, 1e-6] sim.inith = 1e-6 sim.maxh = 1e-6 # Ensure consistent initial conditions sim.make_consistent('IDA_YA_YDP_INIT') # Simulate backward to a small positive time final_time = 0.01 sim.verbosity = 50 t, y, yp = sim.simulate(tfinal=final_time, ncp_list=np.linspace(final_time, bar_omega, 500)) # Retrieve and plot the solution plt.plot(t, y[:, 0], label='Psi') plt.plot(t, y[:, 1], label="Psi'") plt.xlabel('omega') plt.ylabel('Value') plt.legend() plt.title('Backward Solution of HJB DAE') plt.grid(True) plt.show() Continue reading...
