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

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                η→∞  

Answer 1

[Mrutyunjaya Hiremath]

% 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]

for θ and ϕ, initial condition can not be zero.

[Torsten]

Use "bvp4c" or "bvp5c". I don't know which MATLAB code uses the mentionned Keller Box scheme.

[Prachi]

Anyone please give any other suggestions?

[Torsten]

Your question is weird. Is it you who has to implement the Keller Box scheme ?

[Prachi]

I want to apply Keller box scheme for Jeffrey fluid.

[Torsten]

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.

Answer 2

[Santosh Devi]

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

[Santosh Devi]

please give code for this

[Santosh Devi]

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

[Torsten]

If you have code, what is your question ?