I'm trying to solve jeffry hamel equation using RK4 but it's not giving output as expected.

%%%% Jeffry-Hamel equation stability analysis %%%%
%%% ode:(U''' +2sUU' = 0) such that U' = g, g' = h, h' = -2*s*U*g %%%
%%% where s = S/a  and a = slope of the wall %%%
%%% BC: U =1,g =0 at y =0 && U = 0, g= 10 (Kn = 0.1) at y = 1 %%%
%%% Find h(0) such that g(1) = 10 %%%
%% Range of Variable %%  
y_cl = 0;
y_chw = 1;
Re = 0.68e5;
S = 6;
%% Boundary conditions %%
U0 = 1;
U1 = 0;
g0 = 0;
g1 = 10;
h0 = [0.1 0.3];
%% Improved step size and error tolerance %%
d_y = 0.001;  % Reduce step size for better accuracy
N = (y_chw -y_cl)/d_y;
err = 1;
max_iter = 6e5;  % Limit the maximum number of iterations
iter_count = 0;
%% initializing solution %%
y = y_cl : d_y : y_chw;     % y coordinate for plot function
while err > 1e-15 && iter_count < max_iter
    iter_count = iter_count + 1;
    for j = 1:2
        F = zeros(N+1, 3);  % Initialize solution matrix
        F(1, :) = [U0 g0 h0(j)];
        
        % Runge-Kutta 4th order method
        for i = 1:N
            k1 = d_y * [F(i, 2) F(i, 3) -(F(i, 1) * F(i, 2)) * 2 * S];
            k2 = d_y * [F(i, 2) + k1(2)/2 F(i, 3) + k1(3)/2 -( (F(i, 1) + k1(1)/2) * (F(i, 2) + k1(2)/2) ) * 2 * S];
            k3 = d_y * [F(i, 2) + k2(2)/2 F(i, 3) + k2(3)/2 -( (F(i, 1) + k2(1)/2) * (F(i, 2) + k2(2)/2) ) * 2 * S];
            k4 = d_y * [F(i, 2) + k3(2) F(i, 3) + k3(3) -( (F(i, 1) + k3(1)) * (F(i, 2) + k3(2)) ) * 2 * S];
            F(i+1, :) = F(i, :) + (1/6) * (k1 + 2*k2 + 2*k3 + k4);
        end
        
        g_1(j) = F(end, 2);  % Store g(1) for each trial
    end
    
    % Compute error and update h0
    [err, index] = max(abs(g1 - g_1));
    P = diff(g_1);
    if P ~= 0
        h0_new = h0(1) + (diff(h0) / diff(g_1)) * (g1 - g_1(1));
        h0(index) = h0_new;
    end
end
% Final check on iteration count
if iter_count >= max_iter
    warning('Max iterations reached. Convergence may not be achieved.');
end
%% plotting
f1=figure;
f1.Units = 'normalized';
f1.Position = [0.1 0.1 0.8 0.6];
plot(F,y','LineWidth',1);
xlim([0 1]);
ylim([0 1]);
axis square
yticks(0:0.2:1);
yticklabels({'0','0.2','0.4','0.6','0.8','1'});
xticks(0:0.2:1);
xticklabels({'0','0.2','0.4','0.6','0.8','1'});
grid on
title('Jeffry Hamel solution');
xlabel('y');
legend('U','U"','U"''');      

I'm also inserting a plot the I want to see:

please help me in getting this resolved.

Thanks in advance.

Answer 1

[Alan Stevens]

The following code produces something very close to the plot that you want to see. However, the values of g(1) are nowhere near 10. Are you sure you have expressed the Jeffery-Hamel equations correctly?

%%%% Jeffry-Hamel equation stability analysis %%%%
%%% ode:(U''' +2sUU' = 0) such that U' = g, g' = h, h' = -2*s*U*g %%%
%%% where s = S/a  and a = slope of the wall %%%
%%% BC: U =1,g =0 at y =0 && U = 0, g= 10 (Kn = 0.1) at y = 1 %%%
%%% Find h(0) such that g(1) = 10 %%%

yspan = 0:0.1:1;
U0 = 1; g0 = 0; 
h0 = [-1.7, -2.0, -4.3, -4.1];
s = [0, 0, 6, 6];
sl = ['-s'; '-s'; '-+'; '-+'];

figure
hold on

for i = 1:numel(h0)
Ugh0 = [U0, g0, h0(i)];

[y, Ugh] = ode45(@(y,Ugh)fn(y,Ugh,s(i)), yspan, Ugh0);
U = Ugh(:,1); g = Ugh(:,2); h = Ugh(:,3);

plot(U,y,sl(i,:))
disp(['with s = ', num2str(s(i)),' and h(0) = ',num2str(h0(i)),...
      ' g(1) = ', num2str(g(end))])
end

with s = 0 and h(0) = -1.7 g(1) = -1.7
with s = 0 and h(0) = -2 g(1) = -2
with s = 6 and h(0) = -4.3 g(1) = -0.78216
with s = 6 and h(0) = -4.1 g(1) = -0.67057

grid
xlabel('U'), ylabel('y')
axis([0,1,0,1])
legend(['s = 0', '  h(0) = ', num2str(h0(1))], ...
       ['s = 0 ','  h(0) = ', num2str(h0(2))], ...
       ['s = 6', '  h(0) = ', num2str(h0(3))], ...
       ['s = 6', '  h(0) = ', num2str(h0(4))])

function dUghdy = fn(~,Ugh,s)
      
        U = Ugh(1); g = Ugh(2); h = Ugh(3);
        dUghdy = [g; h; -2*s*U*g];

end

[Kushagra Saurabh]

Dear Alan, Thank you, so much for kind help. I was really struggling with it but now I think I could figure it out. Thanks once again!

[Kushagra Saurabh]

Dear Alan,

I have a question, if suppose I want to find the solution to any 's' say 52 then how should i know that the following h0 is correct for that? is there any way that h can change automatically with change in my 's' ?

thanks for your kind help.

[Alan Stevens]

You need to know the target value(s) of some parameter. I simply adjusted the values of h(0) "by eye" until the resulting graphs looked roughly like those in the plot you provided. For example:

yspan = 0:0.01:1;
U0 = 1; g0 = 0; 
h0 = [-20, -30] ;
s = [52, 52]; 
sl = ['-s'; '-s'];

figure
hold on

for i = 1:numel(h0)
Ugh0 = [U0, g0, h0(i)];

[y, Ugh] = ode45(@(y,Ugh)fn(y,Ugh,s(i)), yspan, Ugh0);
U = Ugh(:,1); g = Ugh(:,2); h = Ugh(:,3);

plot(U,y,sl(i,:))
disp(['with s = ', num2str(s(i)),' and h(0) = ',num2str(h0(i)),...
      ' g(1) = ', num2str(g(end))])
end

with s = 52 and h(0) = -20 g(1) = -0.84233
with s = 52 and h(0) = -30 g(1) = -3.018

grid
xlabel('U'), ylabel('y')
axis([0,1,0,1])
legend(['s = 0', '  h(0) = ', num2str(h0(1))], ...
       ['s = 0 ','  h(0) = ', num2str(h0(2))]) ...

function dUghdy = fn(~,Ugh,s)
      
        U = Ugh(1); g = Ugh(2); h = Ugh(3);
        dUghdy = [g; h; -2*s*U*g];

end

But I have no idea if the resulting plots are at all meaningful!

[Kushagra Saurabh]

i got it that by adusting the h0 we can get the results but, my question is can we make h0 such that it adjust itself according to any s? because there is no way i could validate the result for s more than 6 (as shown in plot i provided)

[Alan Stevens]

"...can we make h0 such that it adjust itself according to any s?..."

Yes, as long as you have some sort of target value(s) for some relevant parameter.

[Kushagra Saurabh]

hey! Alan,

I was trying using another function that iterates the guess value and give the results , it worked well for blasius but in my case it not givving the reults as expected , what should i change here?

here's my code:

clc; clear all;

% Parameters

s = 6; % Example value for s = S/a, adjust as needed

% Initial guess for the shooting method

F_init = [1, 0, -5]; % [U(0), g(0), h(0)]

tol = 1e-6; % Tolerance for boundary condition

max_iter = 10000; % Maximum number of iterations

% Solve using the shooting method

[x, y, F_final] = shooting(F_init, s, tol, max_iter);

% Plot the results

figure; hold on;

plot(x, y(:,1), 'k-', 'LineWidth', 2); % U

plot(x, y(:,2), 'r-', 'LineWidth', 2); % g

plot(x, y(:,3), 'b-', 'LineWidth', 2); % h

figformat;

%% Supplementary Functions

function [x, y, F_final] = shooting(F_init, s, tol, max_iter)

% Manually adjust the initial guess for F using iteration

F = F_init;

for iter = 1:max_iter

% Solve the ODE-IVP with the current guess F

[x, y] = solve_ode(F, s);

% Evaluate the boundary condition at y = 1: U(1) should be 0

U_end = y(end, 1); % U(1)

error = U_end - 0; % Error in boundary condition

if abs(error) < tol

fprintf('Converged after %d iterations with error %e\n', iter, error);

break;

else

% Adjust the initial guess (basic update rule)

F(3) = F(3) - 0.1 * error; % Update h(0) based on error

end

end

if iter == max_iter

fprintf('Did not converge after %d iterations\n', max_iter);

end

F_final = F;

end

function [x, y] = solve_ode(F, s)

% Solve the ODE-IVP with initial condition F on [0 1] (boundary at y = 1)

[x, y] = ode45(@(x, y) [y(2); y(3); -2*s*y(1)*y(2)], [0 10], F);

end

function figformat

% This function simply formats the figure

xlabel('y');

ylabel('U, g, h');

xlim([0 1]);

xticks(0:0.2:1);

h = legend('\it{U}','\it{g}','\it{h}','Location','NorthEast');

legend boxoff;

set(h,'FontSize',18);

axis square

set(gca, 'box', 'on');

set(gca,'FontWeight','normal','linewidth',1.05);

fontname = 'Arial';

set(0,'defaultaxesfontname',fontname);

set(0,'defaulttextfontname',fontname);

fontsize = 20;

set(0,'defaultaxesfontsize',fontsize);

set(0,'defaulttextfontsize',fontsize);

end

thanks for your help.

[Alan Stevens]

You don't say what results you are expecting, but changing the integration end time to 1 from 10 produces the following:

% Parameters
s = 6; % Example value for s = S/a, adjust as needed
% Initial guess for the shooting method
F_init = [1, 0, -5]; % [U(0), g(0), h(0)]
tol = 1e-6; % Tolerance for boundary condition
max_iter = 1000; % Maximum number of iterations
% Solve using the shooting method
[x, y, F_final] = shooting(F_init, s, tol, max_iter);

Converged after 33 iterations with error -9.684316e-07

% Plot the results
figure; hold on;
plot(x, y(:,1), 'k-', 'LineWidth', 2); % U
plot(x, y(:,2), 'r-', 'LineWidth', 2); % g
plot(x, y(:,3), 'b-', 'LineWidth', 2); % h
figformat;

%% Supplementary Functions
function [x, y, F_final] = shooting(F_init, s, tol, max_iter)
    % Manually adjust the initial guess for F using iteration
    F = F_init;
    
    for iter = 1:max_iter
        % Solve the ODE-IVP with the current guess F
        [x, y] = solve_ode(F, s);
        
        % Evaluate the boundary condition at y = 1: U(1) should be 0
        U_end = y(end, 1); % U(1)
        error = U_end - 0; % Error in boundary condition
        
        if abs(error) < tol
            fprintf('Converged after %d iterations with error %e\n', iter, error);
            break;
        else
            % Adjust the initial guess (basic update rule)
            F(3) = F(3) -  error; % Update h(0) based on error   Converges quicker here if error not multiplied by a small factor
        end
    end
    
    if iter == max_iter
        fprintf('Did not converge after %d iterations\n', max_iter);
    end
    
    F_final = F;
end
function [x, y] = solve_ode(F, s)
    
    % Solve the ODE-IVP with initial condition F on [0 1] (boundary at y = 1)
    [x, y] = ode45(@(x, y) [y(2); y(3); -2*s*y(1)*y(2)], [0 1], F);  %%%% I've reset the end of the integration time to 1 rather than 10
end
function figformat
    % This function simply formats the figure
    xlabel('y');
    ylabel('U, g, h');
    xlim([0 1]);
    xticks(0:0.2:1);
    
    h = legend('\it{U}','\it{g}','\it{h}','Location','NorthEast');
    legend boxoff;
    set(h,'FontSize',18);
    axis square
    set(gca, 'box', 'on');
    set(gca,'FontWeight','normal','linewidth',1.05);
    fontname = 'Arial';
    set(0,'defaultaxesfontname',fontname);
    set(0,'defaulttextfontname',fontname);
    fontsize = 20;
    set(0,'defaultaxesfontsize',fontsize);
    set(0,'defaulttextfontsize',fontsize);
end

[Kushagra Saurabh]

sorry for not being clear earlier but the boundary conditions are actually like this: U =1, g=0 at y=0 and U = 0 at y = 1 that's why i changed it. I know that this the plot I will be getting by using the code above but i want the plot similar to one we achived earlier(as like the paper) and I also want that my shooting step optimises atomatically, for any random S.

Thanks for your help.

[Alan Stevens]

Can you upload a copy of the paper (or a link to it)?

[Kushagra Saurabh]

here's the paper, which I'm refering to:

[Torsten]

I can only find the boundary condition U + K*n*dU/dy = 0 at y = 1. So you want to set K = 0 ?

Answer 2

[Torsten]

s = 6;
Kn = 0.1;
U20 = -1;
U2 = fsolve(@(U2)fun_shoot(U2,s,Kn),U20)

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

U2 = -4.0971

[Y,Z] = fun_odes(U2,s);
plot(Z(:,1),Y)
xlabel('U')
ylabel('y')
grid on

function res = fun_shoot(U2,s,Kn)
  [Y,Z] = fun_odes(U2,s);
  res = Z(end,1) + Kn* Z(end,2);
end

function [Y,Z] = fun_odes(U2,s)
  f = @(y,z)[z(2);z(3);-2*s*z(1)*z(2)];
  z0 = [1;0;U2];
  [Y,Z] = ode45(f,[0 1],z0);
end

[Kushagra Saurabh]

Thank you so much Torsten, it would be helpful if you coud just help me know where I was wrong so that I can learn.

Thanks for your kind help.

[Torsten]

I don't know - I didn't check your code. I just used MATLAB functions ("fsolve" for shooting and "ode45" for integration) and it worked fine.

So what you should learn is to use as many MATLAB functions and to implement as little self-made numerical parts as possible to solve a problem.

E.g. the update

F(3) = F(3) -  error; % Update h(0) based on error   Converges quicker here if error not multiplied by a small factor

doesn't make sense at all in your code. How should the values you get for U at y=1 influence the 2nd derivative of U at y=0 in this way ?