How do I solve a system using rkf45 method with shooing technique?

I tried to solve a ode system using rk method with shooing technique. But I have a lot of errors. How can I solve this issue.

bvp_shooting_rk_2()

Not enough input arguments.

Error in solution>ode_system (line 102)
f_triple_prime = A1 * ((f_prime)^2 - f * f_double_prime) + A2 * M * f_prime + K_p * f_prime - beta * A1 * (f^2 * f_triple_prime - 2 * f * f_prime * f_double_prime);

Error in solution>@(eta,y)ode_system(eta,y,A1,A2,A3,A4,K_p,beta,Pr,k_f,k_hnf,Q,Ec) (line 64)
[eta, y] = ode45(@(eta, y) ode_system(eta, y, A1, A2, A3, A4, K_p, beta, Pr, k_f, k_hnf, Q, Ec), [0 eta_end], y0, options);

Error in odearguments (line 92)
f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 104)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

Error in solution>shooting_method (line 64)
[eta, y] = ode45(@(eta, y) ode_system(eta, y, A1, A2, A3, A4, K_p, beta, Pr, k_f, k_hnf, Q, Ec), [0 eta_end], y0, options);

Error in solution>bvp_shooting_rk_2 (line 37)
[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);

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 = 25; % 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);

% 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)
% Shooting method to adjust xi
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-6);

% Initial conditions
y0 = [0; 1; 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), [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; 1; 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, 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, M)
% Extract variables from y
f = y(1);
f_prime = y(2);
f_double_prime = y(3);
theta = y(4);
theta_prime = y(5);

% Define the system of ODEs
f_triple_prime = A1 * ((f_prime)^2 - f * f_double_prime) + A2 * M * f_prime + K_p * f_prime - beta * A1 * (f^2 * f_triple_prime - 2 * f * f_prime * f_double_prime);
theta_double_prime = k_f * Pr * (A3 * f * theta_prime - Ec * A4 * (f_double_prime)^2 - Q * theta) / k_hnf;

% Return derivatives as a column vector
dydeta = [f_prime; f_double_prime; f_triple_prime; theta_prime; theta_double_prime];
end

[Prakash]

Kindly resolve the following errors

Error using bvp_shooting_rk_2>ode_system (line 100)

Not enough input arguments.

Error in bvp_shooting_rk_2>@(eta,y)ode_system(eta,y,A1,A2,A3,A4,K_p,beta,Pr,k_f,k_hnf,Q,Ec) (line 62)

[eta, y] = ode45(@(eta, y) ode_system(eta, y, A1, A2, A3, A4, K_p, beta, Pr, k_f, k_hnf, Q, Ec), [0 eta_end], y0, options);

Error in odearguments (line 87)

f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 113)

[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Error in bvp_shooting_rk_2>shooting_method (line 62)

[eta, y] = ode45(@(eta, y) ode_system(eta, y, A1, A2, A3, A4, K_p, beta, Pr, k_f, k_hnf, Q, Ec), [0 eta_end], y0, options);

Error in bvp_shooting_rk_2 (line 35)

[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);

[Torsten]

@Walter Roberson already answered this: You forgot to pass "M" to "ode_system".

[Prakash]

Thankyou for your commends. The current issue has been solved, but shooting method is faild to converges. What can I do? I attached my revised code. Kindly check the code and resolve this issue

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; 1; 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; 1; 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
y(2); 
y(3);
A1 * (y(2)^2 - y(1) * y(3)) - A2 * M * y(2) - K_p * y(2) + beta * A1 * (y(1)^2 * y(4) - 2 * y(1) * y(2) * y(3));
y(5); 
k_f * Pr * (A3 * y(1) * y(5) - Ec * A4 * y(3)^2 - Q * y(5)) / k_hnf

% Return derivatives as a column vector
dydeta = [y(1); y(2); y(3); y(4); y(5)];
end 

[Torsten]

Why do you think the adjustments you make to "xi_new" are suited to make y(end,2) and y(end,4) converge to 0 ? I don't see that.

    xi_new(1) = xi_new(1) - bc_error_f_prime * 0.1;
    xi_new(2) = xi_new(2) - bc_error_theta * 0.1;

[Prakash]

I want to check the guesses are must satisfies the boundary conditions. ie., y(end,2) and y(end,4) converge to 0

[Torsten]

You must use Newton's method to update the guesses for f''(0) and theta'(0) sensefully depending on what you get at eta = eta_end for y(end,2) and y(end,4). This is a system of two nonlinear equations in two unknowns.

If you can't or don't want to program this on your own, why don't you use "bvp4c" for your problem ?

But as you can see below, even the sophisticated MATLAB solver doesn't converge for your set of equations. You should check them again to see if you made an error in their implementation.

bvp_shooting_rk_2()

Warning: Unable to meet the tolerance without using more than 2000 mesh points.
The last mesh of 862 points and the solution are available in the output argument.
The maximum residual is 208.196, while requested accuracy is 0.001.

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)
  etamesh = linspace(0,eta_end,100);
  solinit = bvpinit(etamesh,[0; 1; xi_guess(1); S; xi_guess(2)]);
  bvpfcn = @(x,y)[y(2); y(3);A1 * (y(2)^2 - y(1) * y(3)) - A2 * M * y(2) - K_p * y(2) + beta * A1 * (y(1)^2 * y(4) - 2 * y(1) * y(2) * y(3));y(5); k_f * Pr * (A3 * y(1) * y(5) - Ec * A4 * y(3)^2 - Q * y(5)) / k_hnf];
  bcfcn = @(ya,yb)[ya(1);ya(2)-1;yb(3);ya(4)-S;yb(5)];
  sol = bvp4c(bvpfcn, bcfcn, solinit);
  xi = [sol.y(3,1),sol.y(5,1)];
  y = (sol.y).';
  eta = (sol.x).';
end

[Prakash]

Thankyou for your commends. I will be improve the solutions of this problem with rkf45 method

Answer 1

[Umang Pandey]

Hi Prakash,

You can refer to the following File Exchange submission for implementing RKF45 in MATLAB:

https://www.mathworks.com/matlabcentral/fileexchange/73881-runge-kutta-fehlberg-rkf45

Best,

Umang

Answer 2

[Walter Roberson]

[eta, y] = ode45(@(eta, y) ode_system(eta, y, A1, A2, A3, A4, K_p, beta, Pr, k_f, k_hnf, Q, Ec), [0 eta_end], y0, options);

versus

function dydeta = ode_system(eta, y, A1, A2, A3, A4, K_p, beta, Pr, k_f, k_hnf, Q, Ec, M)

The function definition expects to receive M after the other variables, but you are not passing in M.