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

Question

Kushagra Saurabh 2024-9-4，12:04

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2150039-i-m-trying-to-solve-jeffry-hamel-equation-using-rk4-but-it-s-not-giving-output-as-expected

编辑： Torsten 2024-9-13，12:24

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Alan Stevens 2024-9-6，10:26

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2150039-i-m-trying-to-solve-jeffry-hamel-equation-using-rk4-but-it-s-not-giving-output-as-expected#answer_1512139

在 MATLAB Online 中打开

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

11 个评论
显示 9更早的评论隐藏 9更早的评论

Kushagra Saurabh 2024-9-13，6:25

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

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 2024-9-13，7:55

在 MATLAB Online 中打开

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

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 2024-9-13，8:48

J005_2008_Sahu_EuropeanPhysicalJournal-AppliedPhysics-1.pdf

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

Torsten 2024-9-13，9:55

编辑：Torsten 2024-9-13，9:57

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 2024-9-13，10:18

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2150039-i-m-trying-to-solve-jeffry-hamel-equation-using-rk4-but-it-s-not-giving-output-as-expected#answer_1516085

编辑：Torsten 2024-9-13，10:29

在 MATLAB Online 中打开

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

2 个评论
显示无隐藏无

Kushagra Saurabh 2024-9-13，10:43

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 2024-9-13，10:48

编辑：Torsten 2024-9-13，12:24

在 MATLAB Online 中打开

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 ?

请先登录，再进行评论。

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

11 个评论
显示 9更早的评论隐藏 9更早的评论

更多回答（1 个）

2 个评论
显示无隐藏无

另请参阅

类别

标签

Community Treasure Hunt

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

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

11 个评论 显示 9更早的评论隐藏 9更早的评论

更多回答（1 个）

2 个评论 显示 无隐藏 无

另请参阅

类别

标签

Community Treasure Hunt

WeChat

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

11 个评论
显示 9更早的评论隐藏 9更早的评论

2 个评论
显示无隐藏无