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.%%%% 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. %%%% 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. rk4, matlab, channelflow MATLAB Answers — New Questions
