How can I solve the problems "Shooting method failed to converge"?

I solved a bvp problem using shooting thechnique. But "Shooting method failed to converge". How do I solve this issue and How to check the guesses values are correct to satisfies boundary conditions.

eq1 = diff(f(x), x, x, x) - A1*((diff(f(x), x))^2 - f(x)*diff(f(x), x, x)) - A2*M*diff(f(x), x) - Kp*diff(f(x), x) + beta*A1*(f(x)^2*diff(f(x), x, x, x) - 2*f(x)*diff(f(x), x)*diff(f(x), x, x)) == 0;

eq2 = diff(theta(x), x, x) - kf*Pr*(A3*f(x)*diff(theta(x), x) - Ec*A4*diff(f(x), x, x)^2 - Q*theta(x))/khnf == 0;

ics = [f(0) == s, diff(f(x), x, 1) == 1, khnf*diff(theta(x), x)/kf == Bi*(1-theta(0))];

bcs = [f'(5) == 0, theta(5) == 0];

My codes are:

function bvp_shooting_rk_2()

% Define constants and parameters
rho_f = 783; rCp_f = 2090 * rho_f; k_f = 0.145; rb_f = 21e-6; Pr = 6;
rho_s1 = 1800; rCp_s1 = 717 * rho_s1; k_s1 = 5000; rb_s1 = 6.30e7;
rho_s2 = 10500; rCp_s2 = 235 * rho_s2; k_s2 = 429; rb_s2 = 63e-6;
m = 3; phi_s1 = 0.01; phi_s2 = 0.01; M = 0; K_p = 0.5; beta = 0.5;
S = 0.5; % Initial condition for theta at eta = 0
B_i = 0.5; Q = 0.01; Ec = 0.8;

% Assume values for mu_f and sigma_f, as they are not defined in the problem
mu_f = 1; % Dynamic viscosity of the base fluid (assumed)
sigma_f = 1; % Electrical conductivity of the base fluid (assumed)

rho_hnf = (1-phi_s2)*((1-phi_s1)*rho_f + phi_s1*rho_s1) + phi_s2*rho_s2;
rCp_hnf = (1-phi_s2)*((1-phi_s1)*rCp_f + phi_s1*rCp_s1) + phi_s2*rCp_s2;
mu_hnf = mu_f / ((1-phi_s1)^2.5 * (1-phi_s2)^2.5);
k_bf = k_f * (k_s1 + (m-1)*k_f - (m-1)*phi_s1*(k_f - k_s1)) / (k_s1 + (m-1)*k_f - phi_s1*(k_f - k_s1));
k_hnf = k_bf * (k_s2 + (m-1)*k_bf - (m-1)*phi_s2*(k_bf - k_s2)) / (k_s2 + (m-1)*k_bf - phi_s2*(k_bf - k_s2));
rb_hnf = (1-phi_s2)*((1-phi_s1)*rb_f + phi_s1*rb_s1) + phi_s2*rb_s2;
sigma_bf = sigma_f * (1-phi_s1)^2.5;
sigma_hnf = sigma_bf * (1-phi_s2)^2.5;

A1 = (1-phi_s1)^2.5 * (1-phi_s2)^2.5 * ((1-phi_s2)*(1-phi_s1 + phi_s1*rho_s1/rho_f) + phi_s2*rho_s2/rho_f);
A2 = (1-phi_s2) * (1-phi_s1 + phi_s1*rho_s1/rho_f) + phi_s2*rho_s2/rho_f;
A3 = (1-phi_s2) * (1-phi_s1 + phi_s1*rCp_s1/rCp_f) + phi_s2*rCp_s2/rCp_f;
A4 = (1-phi_s1)^2.5 * (1-phi_s2)^2.5;

% Initial guess for the shooting method
xi_guess = [1; 1]; % Initial guesses for f''(0) and theta'(0)
eta_end = 5; % End of the interval
tol = 1e-6; % Tolerance

% Solve the BVP using the shooting method
[xi, y, eta] = shooting_method(xi_guess, eta_end, tol, A1, A2, A3, A4, K_p, beta, Pr, k_f, k_hnf, Q, Ec, S, M);

% Plot the solution
figure;
subplot(2, 1, 1);
plot(eta, y(:, 1), 'b-', 'LineWidth', 2);
xlabel('\eta');
ylabel('f(\eta)');
title('Solution for f(\eta) using Shooting Method with RK');
grid on;

subplot(2, 1, 2);
plot(eta, y(:, 4), 'r-', 'LineWidth', 2);
xlabel('\eta');
ylabel('\theta(\eta)');
title('Solution for \theta(\eta) using Shooting Method with RK');
grid on;
end

function [xi, y, eta] = shooting_method(xi_guess, eta_end, tol, A1, A2, A3, A4, K_p, beta, Pr, k_f, k_hnf, Q, Ec, S, M)
% Shooting method to adjust xi
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-6);

% Initial conditions
y0 = [0; 0; xi_guess(1); S; xi_guess(2)]; % Initial conditions including S for theta at eta = 0

% Solve the ODE system
[eta, y] = ode45(@(eta, y) ode_system(eta, y, A1, A2, A3, A4, K_p, beta, Pr, k_f, k_hnf, Q, Ec, S, M), [0 eta_end], y0, options);

% Check the boundary condition at eta_end
bc_error_f_prime = y(end, 2);  % f'(eta_end) should be 0
bc_error_theta = y(end, 4);    % theta(eta_end) should be 0

% Adjust xi using a simple iteration (could use a more sophisticated method)
xi_new = xi_guess;
iteration_count = 0;
max_iterations = 100;  % Limit the number of iterations to prevent infinite loops
while (abs(bc_error_f_prime) > tol || abs(bc_error_theta) > tol) && iteration_count < max_iterations
    % Simple adjustment step
    xi_new(1) = xi_new(1) - bc_error_f_prime * 0.1;
    xi_new(2) = xi_new(2) - bc_error_theta * 0.1;
    y0 = [0; 0; xi_new(1); S; xi_new(2)];
    [eta, y] = ode45(@(eta, y) ode_system(eta, y, A1, A2, A3, A4, K_p, beta, Pr, k_f, k_hnf, Q, Ec, S, M), [0 eta_end], y0, options);
    
    bc_error_f_prime = y(end, 2);
    bc_error_theta = y(end, 4);
    iteration_count = iteration_count + 1;
end

if iteration_count >= max_iterations
    error('Shooting method failed to converge.');
end

xi = xi_new;
end

function dydeta = ode_system(eta, y, A1, A2, A3, A4, K_p, beta, Pr, k_f, k_hnf, Q, Ec, S, M)
% Extract variables from y
dy = zeros(5, 1);
dy(1) = y(2); 
dy(2) = y(3);
dy(3) = A1 * (y(2)^2 - y(1) * y(3)) + A2 * M * y(2) + K_p * y(2) - beta * A1 * (y(1)^2 * y(3) - 2 * y(1) * y(2) * y(3));
dy(4) = y(5); 
dy(5) = k_f * Pr * (A3 * y(1) * y(5) - Ec * A4 * y(3)^2 - Q * y(4)) / k_hnf;

% Return derivatives as a column vector
dydeta = dy;
end

Answer 1

[Shubham]

Hi Prakash,

The "Shooting method failed to converge" error generally indicates that the initial guesses for the unknown boundary values are not sufficiently close to the correct values. This issue arises not from MATLAB itself, but from the challenge of identifying appropriate initial conditions within the specified number of iterations.

Some potential workarounds can be:

• Improving your initial estimates for boundary values can significantly enhance the shooting method convergence
• Modifying the step size in accordance with initial estimate, a smaller step size may help in achieving convergence
• Increasing the number of iterations in the shooting method
• Adjusting the tolerance level, a higher tolerance level may help in achieving convergence, but at the cost of accuracy.
• implementing a more sophisticated root-finding technique such as the secant method or Newton-Raphson method rather than using simple iterations
• Adding debugging statements to check the values of the residuals and the shooting parameters after each iteration.

These adjustments can potentially resolve the convergence issue you are encountering.

Hope it helps.