# Can anybody help me to code boundary conditions in MATLAB for Keller Box Method?

10 views (last 30 days)
Prachi on 6 Aug 2023
Commented: Torsten on 27 Feb 2024
Can anybody help me to code boundary conditions in MATLAB for Keller Box Method?
f^'=1,f=S,θ^'=-r_1 [1+θ],ϕ^'=-r_2 [1+ϕ] at η=0
f^'=0,f^''=0,θ=0,ϕ=0 as η→∞

Mrutyunjaya Hiremath on 6 Aug 2023
% Define parameters
r_1 = 0.1;
r_2 = 0.2;
S = 2.0;
% Define the differential equations
% y(1) = f, y(2) = f', y(3) = θ, y(4) = ϕ
ode_system = @(eta, y) [y(2); 1; y(3); y(4)];
% Define the boundary conditions at η = 0
initial_conditions = [S, 1, 0, 0];
% Define the boundary conditions at η → ∞
eta_infinity = 100; % Choose a large value
final_conditions = [0, 0, 0, 0];
% Solve the differential equations
[eta, result] = ode45(ode_system, [0, eta_infinity], initial_conditions);
% Extract the solutions
f = result(:, 1);
f_prime = result(:, 2);
theta = result(:, 3);
phi = result(:, 4);
% Plot the solutions
subplot(2, 2, 1);
plot(eta, f);
xlabel('η');
ylabel('f');
title('f vs. η');
subplot(2, 2, 2);
plot(eta, f_prime);
xlabel('η');
ylabel("f'");
title("f' vs. η");
subplot(2, 2, 3);
plot(eta, theta);
xlabel('η');
ylabel('θ');
title('θ vs. η');
subplot(2, 2, 4);
plot(eta, phi);
xlabel('η');
ylabel('ϕ');
title('ϕ vs. η');
Prachi on 7 Aug 2023
I want to apply Keller box scheme for Jeffrey fluid.
Torsten on 7 Aug 2023
It should be clear that we won't program this for you.
If you have a boundary value problem as above, you can use the MATLAB tools "bvp4c" or "bvp5c".
If your problem is an assignment, you will have to start programming it in MATLAB or make a google search whether you find a MATLAB code that fits your needs.

Santosh Devi on 27 Feb 2024
f^''' (η)+ff^'' (η)-(f^' (η))^2+Mf^' (η)-λf(η)=0
■θ^'' (η)+Pr⁡[f(η)θ^' (η)-b/(u_w^2 ) ηθ(η)+Ec(f^'' (η))^2+Q_0 θ(η)]=0
■ϕ^'' (η)+Sc[f(η) ϕ^' (η)-Kϕ(η)]=0
■f(0)=s,f^' (0)=1,θ(0)=1,ϕ(0)=1
■f^' (∞)→0,θ(∞)→0,ϕ(∞)→0
##### 3 CommentsShow 1 older commentHide 1 older comment
Santosh Devi on 27 Feb 2024
function plot_velocity_vs_eta()
% Parameters
lambda = 0.1;
Pr = 0.72;
b = 1;
uw = 1;
Ec = 0.0001;
Q0 = 1;
Sc = 1;
K = 1;
s = 1;
M_values = [0, 0.5, 1];
% Initialize plot
figure;
hold on;
% Iterate over each Prandtl number
for i = 1:length(M_values)
M = M_values(i);
% Define the system of equations
equations = @(eta, y) [
y(2); % f'(eta) = y(2)
y(3); % f''(eta) = y(3)
-y(1)*y(3) + y(2)^2 - M*y(2) + lambda*y(1); % theta'(eta) = y(3)
-Pr * (y(1)*y(4) - b/(uw^2)*eta*y(3) + Ec*y(3)^2 + Q0*y(4)); % theta''(eta)
y(6); % phi'(eta) = y(6)
-Sc * (y(1)*y(6) - K*y(6)) % phi''(eta)
];
% Define the boundary conditions function
bc = @(ya, yb) boundary_conditions(ya, yb, s);
% Solve the equations
eta_span = [0, 10];
initial_conditions = [s; 0; 1; 1; 0; 0]; % Initial guesses for f, f', theta, phi, f'', phi'
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-9);
[eta, solution] = ode45(equations, eta_span, initial_conditions, options);
% Calculate velocity (f'(eta))
velocity = solution(:, 2);
% Plot velocity against eta
plot(eta, velocity, '-', 'LineWidth', 2, 'DisplayName', sprintf('M = %.2f', M));
end
xlabel('\eta');
ylabel('Velocity (f''(\eta))');
title('Velocity Profile for Different M Numbers');
legend('Location', 'best');
grid on;
hold off;
end
function res = boundary_conditions(ya, yb, s)
% Boundary conditions
res = [ya(1) - s; ya(2); ya(3) - 1; ya(5) - 1; ya(4); ya(6); yb(1) - s; yb(2); yb(3) - 1; yb(5) - 1; yb(4); yb(6)];
end
Torsten on 27 Feb 2024
If you have code, what is your question ?