# shooting method to compute the interface shape

I am trying to use a shooting method to compute the shape of liquid-gas interface given by the following differential equation:

$$\frac{d^2 \theta}{ds^2} = \frac{f(\theta)}{h(h + 3\lambda)}$$

with $$\frac{dh}{ds} = \sin\theta$$. Here, $$f(\theta) = -0.67(\sin\theta)^3/(\theta - \sin \theta \cos \theta)$$, $$\theta$$ is an angle to be solved for, $$s$$ is the parameterization variable and $$\lambda$$ is a given constant. Boundary conditions are: $$\theta|_{s = 0} = \theta_q$$, $$\theta|_{s=l} = \pi/2$$, $$h|_{s=0} = 0$$ and $$h|_{s=l} = d/2$$, where $$\theta_q$$, $$l$$ and $$d$$ are given constants.

I have been told to use $$\theta|_{s=l} = \pi/2$$ as a control condition, i.e., to change the value of $$\theta_q$$ until $$\theta|_{s=l}$$ becomes $$\pi/2$$. Values of $$h|_{s=0}$$ and $$h|_{s=l}$$ are, however, fixed.

To solve this using a shooting method, I denote $$Y = [\theta(s) \ u(s) \ h(s)]^T$$, where $$u = \frac{d\theta}{ds}$$ and rewrite the differential equations as

$$Y' = [u \ \ \frac{f(\theta)}{h(h + 3\lambda)} \ \ \sin \theta]^T$$

I use initial guesses of $$\theta = \theta_q$$ and $$u = 0$$ at $$s = 0$$, solve the ODEs using $$\verb|ode45|$$ in MATLAB and iteratively update the value of $$\theta_q$$ using a bisection method. However, I am unable to get a converged solution and I don't understand where I am doing wrong - see code below for reference. As I am quite new to shooting method, I would appreciate any help in getting this work.

% initialize all parameters
% .....
theta_q_lower = 10*pi/180; % Lower bound of possible theta_q values
theta_q_upper = 68*pi/180; % Upper bound of possible theta_q values

% Shooting method loop
while abs(theta_q_upper - theta_q_lower) > tolerance
% Guess of theta_q
theta_q_guess = (theta_q_lower + theta_q_upper) / 2;

% Solve the ODEs with the current guess for theta_q
[s, Y] = ode45(@(s, Y) odefunc(s, Y, Ca, c, lambda, M), [0, hf], [theta_q_guess; 0; h0]);

% Get the final value of theta from the ODE solution
theta_l = Y(end, 1);

% Adjust the guess based on the result
if theta_l < theta_l_target % theta_q was too low
theta_q_lower = theta_q_guess;
else % theta_q was too high
theta_q_upper = theta_q_guess;
end
iter = iter+1;
end

% Use the final guess for theta_q
theta_q_final = (theta_q_lower + theta_q_upper) / 2;

% Solve the ODEs one last time with the final theta_q
[s, Y] = ode45(@(s, Y) odefunc(s, Y, Ca, c, lambda, M), [0, hf], [theta_q_final; 0; h0]);

% Plot the results
plot(s, Y(:, 1)); % Plot theta(s)
hold on;
plot(s, Y(:, 3)); % Plot h(s)
legend('\theta(s)', 'h(s)');
xlabel('s');
ylabel('Value');
title('Shooting Method Solution');
hold off;

function dYds = odefunc(s, Y, Ca, c, lambda, M)
% Unpack current values
theta = Y(1);
u = Y(2);
h = Y(3);

% Prevent division by zero by ensuring h is not too close to zero
epsilon = 1e-7;
h_safe = h + epsilon;

% Equations
dtheta_ds = u;
du_ds = F(theta, M) / (h_safe * (h_safe + c * lambda));
dh_ds = sin(theta);

dYds = [dtheta_ds; du_ds; dh_ds];
end

function Fval = F(theta, M)
Fval = -0.67 * (sin(theta)^3) / (theta - sin(theta) * cos(theta));
end

Given your differential equation is nonlinear, it is not clear that using the bisection method to find the desired value for $$\theta_q$$ so that $$\theta|_{s = l} = \pi/2$$ is reasonable.
Let $$f(\theta_q)$$ be the value of $$\theta|_{s = l}$$ obtained after using $$\theta|_{s = 0} = \theta_q$$ as the initial condition for solving the differential equation, all other things being fixed. Let us also assume there actually does exist some value $$\theta^*$$ such that $$f(\theta^*) = \pi/2$$. If $$f(\theta_q)$$ is monotonic with respect to $$\theta_q$$, then your bisection approach should work given you do not have issues in your implementation. However, if $$f(\theta_q)$$ is not monotonic, things become trickier. In the worst case that $$f(\theta_q)$$ is nonconvex, then your approach may very well fail to find the desired $$\theta^*$$.
To help diagnose your situation, it may help to evaluate $$f(\theta_q)$$ for many values of $$\theta_q$$ in the range $$[\pi/18, 68 \pi/ 180]$$ and plot a figure to see what it looks like.
• thanks for the suggestion. plot of $f(\theta_q)$ is unfortunately monotonic. Jan 18 at 22:51
• @SthavishthaBhopalam given this is the case, then your binary search algorithm must be wrong because a correct implementation would return the desired answer. If you look at your algorithm, it assumes $f(\theta_q)$ is monotonically increasing. Is this a correct property of $f(\theta_q)$? If it is monotonically decreasing, then your algo will not return a correct result. Arguably, you should write your algorithm so it can handle both cases. Also, when you want the final guess for $\theta_q$, you should use a linear interpolation instead of just taking the average of the final bounds. Jan 20 at 19:29