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;
    iter = iter+1;

% 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)');
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];

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

1 Answer 1


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.

  • $\begingroup$ thanks for the suggestion. plot of $f(\theta_q)$ is unfortunately monotonic. $\endgroup$ Commented Jan 18 at 22:51
  • $\begingroup$ @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. $\endgroup$
    – spektr
    Commented Jan 20 at 19:29

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.