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

10 views (last 30 days)
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 η→∞

Answers (2)

Mrutyunjaya Hiremath
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. η');
  6 Comments
Torsten
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.

Sign in to comment.


Santosh Devi
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 Comments
Santosh Devi
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
% Add labels and legend
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

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!