@contradict The issue with this problem lies in some divergences with the initial values, which require an expansion to resolve. I initially left out the expansion for simplicity, but I’ll add it now.
Regarding the MethodError, this was caused by using an older version of Julia. After updating to 1.10, the error disappeared.
@ChrisRackauckas I’ve been using the least squares method, but it’s highly sensitive to the initial conditions of the free parameters, often leading to different local minima. This is one of the reasons I’m exploring an alternative approach using BVProblem!.
The current issue with the BVProblem approach is that I’m encountering a rank deficiency matrix error. The BVP seems to interpret the problem as having 3 boundary conditions but 6 unknowns (4 initial conditions + 2 free parameters). In reality, only one initial condition, u0, is unknown, making the problem well-posed with 3 boundary conditions and 3 unknowns (1 initial condition and 2 free parameters).
My questions are:
- How can I specify to
BVProblem that not all initial conditions are unknown, and that three of them should remain constant (or be evaluated at each integration call, as I did with the least squares method in the residual function)?
- If this isn’t possible, would it be more convenient to stick with the least squares method? If so, do you have any suggestions for managing the instability of the parameters found?
Using Least squares method:
using DifferentialEquations
using Optim
using Interpolations
function ShapeIntegrator!(dy, y, p, t)
omega,sigma = p
k = 1.0
# dy[1] = DPsi_DS(omega, u)
dy[1] = omega * y[2]
a1 = -y[2] / y[4] * cos(y[1])
a2 = (sin(y[1]) * cos(y[1])) / y[4]^2
a3 = (sin(y[1]) * y[3]) / (2π * y[4] * k)
dy[2] = omega * (a1 + a2 + a3)
b1 = π * k * (y[2]^2 - sin(y[1])^2 / y[4]^2)
b2 = 2π * sigma * k
dy[3] = omega * (b1 + b2)
dy[4] = omega * cos(y[1])
end
# South Pole expansion
function South_X(s, omega, u0)
x1 = omega
x3 = -omega^3 * u0^2
return x1 * s + 1/6 * x3 * s^3
end
function South_Psi(s, omega, u0, sigma)
psi1 = omega * u0
psi3 = omega^3 * (3 * sigma * u0 - 2 * u0^3)
return psi1 * s + 1/6 * psi3 * s^3
end
function South_U(s, omega, u0, sigma)
return u0 + 1/2 * omega^2 * 0.5 * (3 * u0 * sigma - 2 * u0^3) * s^2
end
function South_Gamma(s, omega, u0, sigma, k)
gamma1 = omega * (2π * sigma * k)
gamma3 = 4/3 * π * k * u0 * omega^3 * (3 * sigma * u0 - 2 * u0^3)
return gamma1 * s + 1/6 * gamma3 * s^3
end
function InitialArcLength(p)
# Find the minimimum length s_init to not have a divergence in x(s).
omega = p[1]
u0 = p[2]
threshold = 0.035 # Changed from 0.01 for better stability
n = 1
delta_s = 0.0001
while true
if South_X(n * delta_s, omega, u0) > threshold
break
end
n += 1
end
return n * delta_s
end
function InitialValues(s_init, p)
omega = p[1]
u0 = p[2]
sigma = p[3]
k = p[4]
return [South_Psi(s_init, omega, u0, sigma),
South_U(s_init, omega, u0, sigma),
South_Gamma(s_init, omega, u0, sigma, k),
South_X(s_init, omega, u0),
]
end
function Residuales(parameters, boundary_conditions)
k = 1
omega, sigma, u0 = parameters
# Calculate initial arc length
s_init = InitialArcLength((omega, u0))
# Calculate initial values
z_init = InitialValues(s_init, (omega, u0, sigma, k))
# Evaluation points for the solution
s = range(s_init, stop=1, length=1000)
# Solve the IVP using Radau method
f = ODEFunction(ShapeIntegrator!,jac=ShapeJacobian!)
prob = ODEProblem(f, z_init, (s_init, 1), (omega, sigma))
sol = solve(prob, saveat=s)
if sol.retcode != :Success
println("integration failed")
return [1e3, 1e3, 1e3]
end
z_fina_num = sol.u[end]
psif, uf, xf = boundary_conditions
psi = z_fina_num[1] - psif
u = z_fina_num[2] - uf
x = z_fina_num[4] - xf
return [psi, u, x]
end
function main()
# Constitutive relations
Rparticle = 3
Rvesicle = 30
rpa = Rparticle / Rvesicle
# Wrapping angle
deg = 60
phi = deg2rad(deg)
# Free parameters initial values
omega = 3.0
sigma = 0.019
u0 = 1.0
# Boundary conditions
psistar = π + phi
ustar = 1 / rpa
xstar = rpa * sin(phi)
# Check on xstar value
if xstar < 0.035
println("xstar: $xstar")
error("xstar too small, can lead to divergences")
else
println("xstar: $xstar")
end
boundary_conditions = [psistar, ustar, xstar]
free_params_extended = [omega, sigma, u0]
# Shooting algorithm and solver
result = optimize(x -> sum(abs2, Residuales(x, boundary_conditions)), free_params_extended, NelderMead())
@printf "Err: %.5f\n" result.minimum
best_parameters = result.minimizer
print(best_parameters)
Using BVProblem approach:
function ShapeIntegrator!(dy, y, p, t)
# psistar, ustar, xstar = p
omega = y[5]
sigma = y[6]
dy[1] = omega * y[2]
a1 = -y[2] / y[4] * cos(y[1])
a2 = (sin(y[1]) * cos(y[1])) / y[4]^2
a3 = (sin(y[1]) * y[3]) / (2π * y[4])
dy[2] = omega * (a1 + a2 + a3)
b1 = π * (y[2]^2 - sin(y[1])^2 / y[4]^2)
b2 = 2π * sigma
dy[3] = omega * (b1 + b2)
dy[4] = omega * cos(y[1])
dy[5] = 0 #omega
dy[6] = 0 #sigma
end
function bc1!(residual, y, p, t)
psistar, ustar, xstar = p
residual[1] = y[end][1] - psistar
residual[2] = y[end][2] - ustar
residual[4] = y[end][4] - xstar
end
# South Pole expansion
function South_X(s, omega, u0)
x1 = omega
x3 = -omega^3 * u0^2
return x1 * s + 1/6 * x3 * s^3
end
function South_Psi(s, omega, u0, sigma)
psi1 = omega * u0
psi3 = omega^3 * (3 * sigma * u0 - 2 * u0^3)
return psi1 * s + 1/6 * psi3 * s^3
end
function South_U(s, omega, u0, sigma)
return u0 + 1/2 * omega^2 * 0.5 * (3 * u0 * sigma - 2 * u0^3) * s^2
end
function South_Gamma(s, omega, u0, sigma, k)
gamma1 = omega * (2π * sigma * k)
gamma3 = 4/3 * π * k * u0 * omega^3 * (3 * sigma * u0 - 2 * u0^3)
return gamma1 * s + 1/6 * gamma3 * s^3
end
function InitialArcLength(p)
# Find the minimimum length s_init to not have a divergence in x(s).
omega = p[1]
u0 = p[2]
threshold = 0.035 # Changed from 0.01 for better stability
n = 1
delta_s = 0.0001
while true
if South_X(n * delta_s, omega, u0) > threshold
break
end
n += 1
end
return n * delta_s
end
function InitialValues(s_init, p)
omega = p[1]
u0 = p[2]
sigma = p[3]
k = p[4]
return [South_Psi(s_init, omega, u0, sigma),
South_U(s_init, omega, u0, sigma),
South_Gamma(s_init, omega, u0, sigma, k),
South_X(s_init, omega, u0),
]
end
function test()
# constitutive relations
Rparticle = 3
Rvesicle = 30
rpa = Rparticle / Rvesicle
# wrapping angle
phi = pi/6
# Boundary conditions
psistar = pi + phi
ustar = 1/rpa
xstar = rpa*sin(phi)
p = [psistar, ustar, xstar]
# check on xstar value
if xstar < 0.035
print("xstar:",xstar)
throw(DomainError())
else
print("xstar:",xstar)
end
# free params
omega0 = 3.0
sigma0 = 0.019
u0 = 1.0
k = 1
s_init = InitialArcLength((omega0, u0))
s_span = (s_init, 1.0)
# init conditions
init = InitialValues(s_init, (omega0, u0, sigma0, k))
init = [init; [omega0, sigma0]]
bvp2 = BVProblem(ShapeIntegrator!, bc1!, init, s_span, p)
end
bvp2 = test()
soln = solve(bvp2)