How do I solve a system using rkf45 method with shooing technique?

10 visualizaciones (últimos 30 días)
Prakash
Prakash el 22 de Jun. de 2024
Comentada: Prakash el 26 de Jun. de 2024
I tried to solve a ode system using rk method with shooing technique. But I have a lot of errors. How can I solve this issue.
bvp_shooting_rk_2()
Not enough input arguments.

Error in solution>ode_system (line 102)
f_triple_prime = A1 * ((f_prime)^2 - f * f_double_prime) + A2 * M * f_prime + K_p * f_prime - beta * A1 * (f^2 * f_triple_prime - 2 * f * f_prime * f_double_prime);

Error in solution>@(eta,y)ode_system(eta,y,A1,A2,A3,A4,K_p,beta,Pr,k_f,k_hnf,Q,Ec) (line 64)
[eta, y] = ode45(@(eta, y) ode_system(eta, y, A1, A2, A3, A4, K_p, beta, Pr, k_f, k_hnf, Q, Ec), [0 eta_end], y0, options);

Error in odearguments (line 92)
f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 104)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

Error in solution>shooting_method (line 64)
[eta, y] = ode45(@(eta, y) ode_system(eta, y, A1, A2, A3, A4, K_p, beta, Pr, k_f, k_hnf, Q, Ec), [0 eta_end], y0, options);

Error in solution>bvp_shooting_rk_2 (line 37)
[xi, y, eta] = shooting_method(xi_guess, eta_end, tol, A1, A2, A3, A4, K_p, beta, Pr, k_f, k_hnf, Q, Ec, S);
function bvp_shooting_rk_2()
% Define constants and parameters
rho_f = 783; rCp_f = 2090 * rho_f; k_f = 0.145; rb_f = 21e-6; Pr = 6;
rho_s1 = 1800; rCp_s1 = 717 * rho_s1; k_s1 = 5000; rb_s1 = 6.30e7;
rho_s2 = 10500; rCp_s2 = 235 * rho_s2; k_s2 = 429; rb_s2 = 63e-6;
m = 3; phi_s1 = 0.01; phi_s2 = 0.01; M = 0; K_p = 0.5; beta = 0.5;
S = 0.5; % Initial condition for theta at eta = 0
B_i = 0.5; Q = 0.01; Ec = 0.8;
% Assume values for mu_f and sigma_f, as they are not defined in the problem
mu_f = 1; % Dynamic viscosity of the base fluid (assumed)
sigma_f = 1; % Electrical conductivity of the base fluid (assumed)
rho_hnf = (1-phi_s2)*((1-phi_s1)*rho_f + phi_s1*rho_s1) + phi_s2*rho_s2;
rCp_hnf = (1-phi_s2)*((1-phi_s1)*rCp_f + phi_s1*rCp_s1) + phi_s2*rCp_s2;
mu_hnf = mu_f / ((1-phi_s1)^2.5 * (1-phi_s2)^2.5);
k_bf = k_f * (k_s1 + (m-1)*k_f - (m-1)*phi_s1*(k_f - k_s1)) / (k_s1 + (m-1)*k_f - phi_s1*(k_f - k_s1));
k_hnf = k_bf * (k_s2 + (m-1)*k_bf - (m-1)*phi_s2*(k_bf - k_s2)) / (k_s2 + (m-1)*k_bf - phi_s2*(k_bf - k_s2));
rb_hnf = (1-phi_s2)*((1-phi_s1)*rb_f + phi_s1*rb_s1) + phi_s2*rb_s2;
sigma_bf = sigma_f * (1-phi_s1)^2.5;
sigma_hnf = sigma_bf * (1-phi_s2)^2.5;
A1 = (1-phi_s1)^2.5 * (1-phi_s2)^2.5 * ((1-phi_s2)*(1-phi_s1 + phi_s1*rho_s1/rho_f) + phi_s2*rho_s2/rho_f);
A2 = (1-phi_s2) * (1-phi_s1 + phi_s1*rho_s1/rho_f) + phi_s2*rho_s2/rho_f;
A3 = (1-phi_s2) * (1-phi_s1 + phi_s1*rCp_s1/rCp_f) + phi_s2*rCp_s2/rCp_f;
A4 = (1-phi_s1)^2.5 * (1-phi_s2)^2.5;
% Initial guess for the shooting method
xi_guess = [1; 1]; % Initial guesses for f''(0) and theta'(0)
eta_end = 25; % End of the interval
tol = 1e-6; % Tolerance
% Solve the BVP using the shooting method
[xi, y, eta] = shooting_method(xi_guess, eta_end, tol, A1, A2, A3, A4, K_p, beta, Pr, k_f, k_hnf, Q, Ec, S);
% Plot the solution
figure;
subplot(2, 1, 1);
plot(eta, y(:, 1), 'b-', 'LineWidth', 2);
xlabel('\eta');
ylabel('f(\eta)');
title('Solution for f(\eta) using Shooting Method with RK');
grid on;
subplot(2, 1, 2);
plot(eta, y(:, 4), 'r-', 'LineWidth', 2);
xlabel('\eta');
ylabel('\theta(\eta)');
title('Solution for \theta(\eta) using Shooting Method with RK');
grid on;
end
function [xi, y, eta] = shooting_method(xi_guess, eta_end, tol, A1, A2, A3, A4, K_p, beta, Pr, k_f, k_hnf, Q, Ec, S)
% Shooting method to adjust xi
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-6);
% Initial conditions
y0 = [0; 1; xi_guess(1); S; xi_guess(2)]; % Initial conditions including S for theta at eta = 0
% Solve the ODE system
[eta, y] = ode45(@(eta, y) ode_system(eta, y, A1, A2, A3, A4, K_p, beta, Pr, k_f, k_hnf, Q, Ec), [0 eta_end], y0, options);
% Check the boundary condition at eta_end
bc_error_f_prime = y(end, 2); % f'(eta_end) should be 0
bc_error_theta = y(end, 4); % theta(eta_end) should be 0
% Adjust xi using a simple iteration (could use a more sophisticated method)
xi_new = xi_guess;
iteration_count = 0;
max_iterations = 100; % Limit the number of iterations to prevent infinite loops
while (abs(bc_error_f_prime) > tol || abs(bc_error_theta) > tol) && iteration_count < max_iterations
% Simple adjustment step
xi_new(1) = xi_new(1) - bc_error_f_prime * 0.1;
xi_new(2) = xi_new(2) - bc_error_theta * 0.1;
y0 = [0; 1; xi_new(1); S; xi_new(2)];
[eta, y] = ode45(@(eta, y) ode_system(eta, y, A1, A2, A3, A4, K_p, beta, Pr, k_f, k_hnf, Q, Ec, M), [0 eta_end], y0, options);
bc_error_f_prime = y(end, 2);
bc_error_theta = y(end, 4);
iteration_count = iteration_count + 1;
end
if iteration_count >= max_iterations
error('Shooting method failed to converge.');
end
xi = xi_new;
end
function dydeta = ode_system(eta, y, A1, A2, A3, A4, K_p, beta, Pr, k_f, k_hnf, Q, Ec, M)
% Extract variables from y
f = y(1);
f_prime = y(2);
f_double_prime = y(3);
theta = y(4);
theta_prime = y(5);
% Define the system of ODEs
f_triple_prime = A1 * ((f_prime)^2 - f * f_double_prime) + A2 * M * f_prime + K_p * f_prime - beta * A1 * (f^2 * f_triple_prime - 2 * f * f_prime * f_double_prime);
theta_double_prime = k_f * Pr * (A3 * f * theta_prime - Ec * A4 * (f_double_prime)^2 - Q * theta) / k_hnf;
% Return derivatives as a column vector
dydeta = [f_prime; f_double_prime; f_triple_prime; theta_prime; theta_double_prime];
end
  7 comentarios
Torsten
Torsten el 25 de Jun. de 2024
Editada: Torsten el 25 de Jun. de 2024
You must use Newton's method to update the guesses for f''(0) and theta'(0) sensefully depending on what you get at eta = eta_end for y(end,2) and y(end,4). This is a system of two nonlinear equations in two unknowns.
If you can't or don't want to program this on your own, why don't you use "bvp4c" for your problem ?
But as you can see below, even the sophisticated MATLAB solver doesn't converge for your set of equations. You should check them again to see if you made an error in their implementation.
bvp_shooting_rk_2()
Warning: Unable to meet the tolerance without using more than 2000 mesh points.
The last mesh of 862 points and the solution are available in the output argument.
The maximum residual is 208.196, while requested accuracy is 0.001.
function bvp_shooting_rk_2()
% Define constants and parameters
rho_f = 783; rCp_f = 2090 * rho_f; k_f = 0.145; rb_f = 21e-6; Pr = 6;
rho_s1 = 1800; rCp_s1 = 717 * rho_s1; k_s1 = 5000; rb_s1 = 6.30e7;
rho_s2 = 10500; rCp_s2 = 235 * rho_s2; k_s2 = 429; rb_s2 = 63e-6;
m = 3; phi_s1 = 0.01; phi_s2 = 0.01; M = 0; K_p = 0.5; beta = 0.5;
S = 0.5; % Initial condition for theta at eta = 0
B_i = 0.5; Q = 0.01; Ec = 0.8;
% Assume values for mu_f and sigma_f, as they are not defined in the problem
mu_f = 1; % Dynamic viscosity of the base fluid (assumed)
sigma_f = 1; % Electrical conductivity of the base fluid (assumed)
rho_hnf = (1-phi_s2)*((1-phi_s1)*rho_f + phi_s1*rho_s1) + phi_s2*rho_s2;
rCp_hnf = (1-phi_s2)*((1-phi_s1)*rCp_f + phi_s1*rCp_s1) + phi_s2*rCp_s2;
mu_hnf = mu_f / ((1-phi_s1)^2.5 * (1-phi_s2)^2.5);
k_bf = k_f * (k_s1 + (m-1)*k_f - (m-1)*phi_s1*(k_f - k_s1)) / (k_s1 + (m-1)*k_f - phi_s1*(k_f - k_s1));
k_hnf = k_bf * (k_s2 + (m-1)*k_bf - (m-1)*phi_s2*(k_bf - k_s2)) / (k_s2 + (m-1)*k_bf - phi_s2*(k_bf - k_s2));
rb_hnf = (1-phi_s2)*((1-phi_s1)*rb_f + phi_s1*rb_s1) + phi_s2*rb_s2;
sigma_bf = sigma_f * (1-phi_s1)^2.5;
sigma_hnf = sigma_bf * (1-phi_s2)^2.5;
A1 = (1-phi_s1)^2.5 * (1-phi_s2)^2.5 * ((1-phi_s2)*(1-phi_s1 + phi_s1*rho_s1/rho_f) + phi_s2*rho_s2/rho_f);
A2 = (1-phi_s2) * (1-phi_s1 + phi_s1*rho_s1/rho_f) + phi_s2*rho_s2/rho_f;
A3 = (1-phi_s2) * (1-phi_s1 + phi_s1*rCp_s1/rCp_f) + phi_s2*rCp_s2/rCp_f;
A4 = (1-phi_s1)^2.5 * (1-phi_s2)^2.5;
% Initial guess for the shooting method
xi_guess = [1; 1]; % Initial guesses for f''(0) and theta'(0)
eta_end = 5; % End of the interval
tol = 1e-6; % Tolerance
% Solve the BVP using the shooting method
[xi, y, eta] = shooting_method(xi_guess, eta_end, tol, A1, A2, A3, A4, K_p, beta, Pr, k_f, k_hnf, Q, Ec, S, M);
% Plot the solution
figure;
subplot(2, 1, 1);
plot(eta, y(:,1), 'b-', 'LineWidth', 2);
xlabel('\eta');
ylabel('f(\eta)');
title('Solution for f(\eta) using Shooting Method with RK');
grid on;
subplot(2, 1, 2);
plot(eta, y(:,4), 'r-', 'LineWidth', 2);
xlabel('\eta');
ylabel('\theta(\eta)');
title('Solution for \theta(\eta) using Shooting Method with RK');
grid on;
end
function [xi, y, eta] = shooting_method(xi_guess, eta_end, tol, A1, A2, A3, A4, K_p, beta, Pr, k_f, k_hnf, Q, Ec, S, M)
etamesh = linspace(0,eta_end,100);
solinit = bvpinit(etamesh,[0; 1; xi_guess(1); S; xi_guess(2)]);
bvpfcn = @(x,y)[y(2); y(3);A1 * (y(2)^2 - y(1) * y(3)) - A2 * M * y(2) - K_p * y(2) + beta * A1 * (y(1)^2 * y(4) - 2 * y(1) * y(2) * y(3));y(5); k_f * Pr * (A3 * y(1) * y(5) - Ec * A4 * y(3)^2 - Q * y(5)) / k_hnf];
bcfcn = @(ya,yb)[ya(1);ya(2)-1;yb(3);ya(4)-S;yb(5)];
sol = bvp4c(bvpfcn, bcfcn, solinit);
xi = [sol.y(3,1),sol.y(5,1)];
y = (sol.y).';
eta = (sol.x).';
end
Prakash
Prakash el 26 de Jun. de 2024
Thankyou for your commends. I will be improve the solutions of this problem with rkf45 method

Iniciar sesión para comentar.

Respuestas (2)

Umang Pandey
Umang Pandey el 23 de Jun. de 2024
Hi Prakash,
You can refer to the following File Exchange submission for implementing RKF45 in MATLAB:
Best,
Umang

Walter Roberson
Walter Roberson el 24 de Jun. de 2024
[eta, y] = ode45(@(eta, y) ode_system(eta, y, A1, A2, A3, A4, K_p, beta, Pr, k_f, k_hnf, Q, Ec), [0 eta_end], y0, options);
versus
function dydeta = ode_system(eta, y, A1, A2, A3, A4, K_p, beta, Pr, k_f, k_hnf, Q, Ec, M)
The function definition expects to receive M after the other variables, but you are not passing in M.

Categorías

Más información sobre Marine and Underwater Vehicles en Help Center y File Exchange.

Productos


Versión

R2015a

Community Treasure Hunt

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

Start Hunting!

Translated by