Unable to meet integration tolerances
3 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Hello guys , I can't seem to find a solution on how to go about this warning . Any suggestions are highly appreciated .
Attached below is my code .
close all
clear
clc
% Constants
Pr_air = 0.7;
Pr_water = 6.7;
eta_end = 10; % Assuming this value is sufficiently large to approximate +infinity
delta_eta = 0.05;
% Initial conditions
y0_air = [0, 0, 0.33206, 1, 0]; % f(0), f'(0), f''(0), theta(0), theta'(0) for air
y0_water = [0, 0, 0.33206, 1, 0]; % f(0), f'(0), f''(0), theta(0), theta'(0) for water
% Solve for air and water
sol_air = solve_with_adjustment(Pr_air, y0_air, eta_end, delta_eta);
sol_water = solve_with_adjustment(Pr_water, y0_water, eta_end, delta_eta);
% Plot the results
figure;
% Velocity profile
subplot(1, 2, 1);
plot(sol_air.eta, sol_air.y(:, 2), 'DisplayName', sprintf('Air (Pr=%.1f)', Pr_air));
hold on;
plot(sol_water.eta, sol_water.y(:, 2), 'DisplayName', sprintf('Water (Pr=%.1f)', Pr_water));
title('Similarity Velocity Distribution f''(\eta)');
xlabel('\eta');
ylabel('f''(\eta)');
legend;
hold off;
% Temperature profile
subplot(1, 2, 2);
plot(sol_air.eta, sol_air.y(:, 4), 'DisplayName', sprintf('Air (Pr=%.1f)', Pr_air));
hold on;
plot(sol_water.eta, sol_water.y(:, 4), 'DisplayName', sprintf('Water (Pr=%.1f)', Pr_water));
title('Similarity Temperature Distribution \theta(\eta)');
xlabel('\eta');
ylabel('\theta(\eta)');
legend;
hold off;
% Local functions
function sol = solve_with_adjustment(Pr, y0, eta_end, delta_eta)
% Function to execute the Runge-Kutta integration and adjust f''(0) iteratively
% Initialize variables for the iteration
f_double_prime_0_lower = 0; % Lower bound for f''(0)
f_double_prime_0_upper = 1; % Upper bound for f''(0)
tolerance = 1e-4; % Tolerance for f'(inf)
max_iterations = 100; % Maximum number of iterations
options = odeset('MaxStep', delta_eta);
for i = 1:max_iterations
% Solve the ODEs
y0(3) = (f_double_prime_0_lower + f_double_prime_0_upper) / 2; % Update the guess for f''(0)
[eta_vals, y_vals] = ode45(@(eta, y) odes_system(eta, y, Pr), [0, eta_end], y0, options);
% Check the boundary condition at infinity
if abs(y_vals(end, 2)) < tolerance % If f'(eta_end) is close enough to 0
sol.eta = eta_vals;
sol.y = y_vals;
return;
elseif y_vals(end, 2) > 0 % f'(eta_end) is too high, decrease f''(0)
f_double_prime_0_upper = y0(3);
else % f'(eta_end) is too low, increase f''(0)
f_double_prime_0_lower = y0(3);
end
end
error('Solution did not converge within the maximum number of iterations');
end
function dy = odes_system(eta, y, Pr)
% System of ODEs
dy = zeros(5,1);
dy(1) = y(2);
dy(2) = y(3);
dy(3) = (4/5)*y(2)^2 - (12/5)*y(1)*y(3) - y(4);
dy(4) = y(5);
dy(5) = -(12/5)*Pr*(y(2)*y(4) + y(1)*y(5));
end
2 commentaires
Réponse acceptée
Torsten
le 4 Fév 2024
Modifié(e) : Torsten
le 4 Fév 2024
I reduced Pr_water by a factor of 0.6 to achieve convergence.
% Constants
Pr_air = 0.7;
Pr_water = 6.7*0.6;
eta_end = 4; % Assuming this value is sufficiently large to approximate +infinity
delta_eta = 0.05;
% Initial conditions
y0_air = [0, 0, 0.33206, 1, 0]; % f(0), f'(0), f''(0), theta(0), theta'(0) for air
y0_water = [0, 0, 0.33206, 1, 0]; % f(0), f'(0), f''(0), theta(0), theta'(0) for water
% Solve for air
bvpfcn = @(eta,y)[y(2);y(3);(4/5)*y(2)^2-(12/5)*y(1)*y(3)-y(4);y(5);-(12/5)*Pr_air*(y(2)*y(4)+y(1)*y(5))];
bcfcn = @(etaa,etab)[etaa(1);etaa(2);etaa(4)-1;etaa(5);etab(2)];
etamesh = linspace(0,eta_end,1000);
solinit = bvpinit(etamesh, [0; 0; 0.33206; 1; 0]);
sol_air = bvp4c(bvpfcn, bcfcn, solinit);
% Solve for water
bvpfcn = @(eta,y)[y(2);y(3);(4/5)*y(2)^2-(12/5)*y(1)*y(3)-y(4);y(5);-(12/5)*Pr_water*(y(2)*y(4)+y(1)*y(5))];
sol_water = bvp4c(bvpfcn, bcfcn, solinit);
figure(1)
plot(sol_air.x, sol_air.y(2,:), '-o')
hold on
plot(sol_water.x, sol_water.y(2,:), '-o')
hold off
grid on
figure(2)
plot(sol_air.x, sol_air.y(4,:), '-o')
hold on
plot(sol_water.x, sol_water.y(4,:), '-o')
hold off
grid on
3 commentaires
Torsten
le 4 Fév 2024
Modifié(e) : Torsten
le 4 Fév 2024
I still use Pr as a constant, but had to reduce its value for water from 6.7 to 0.6*6.7 to achieve convergence of the boundary value solver.
Maybe there is a MATLAB version of the python ODE solver RK45 that you could use for your code:
A gradual increase of Pr_water seems to solve the issue:
% Constants
Pr_air = 0.7;
Pr_water = 6.7;
eta_end = 4; % Assuming this value is sufficiently large to approximate +infinity
% Initial conditions
y0_air = [0, 0, 0.33206, 1, 0]; % f(0), f'(0), f''(0), theta(0), theta'(0) for air
y0_water = [0, 0, 0.33206, 1, 0]; % f(0), f'(0), f''(0), theta(0), theta'(0) for water
% Solve for air
bvpfcn = @(eta,y)[y(2);y(3);(4/5)*y(2)^2-(12/5)*y(1)*y(3)-y(4);y(5);-(12/5)*Pr_air*(y(2)*y(4)+y(1)*y(5))];
bcfcn = @(etaa,etab)[etaa(1);etaa(2);etaa(4)-1;etaa(5);etab(2)];
etamesh = linspace(0,eta_end,1000);
solinit = bvpinit(etamesh, [0; 0; 0.33206; 1; 0]);
sol_air = bvp4c(bvpfcn, bcfcn, solinit);
% Solve for water
for i = 1:5
Pr = Pr_water*(0.6+0.1*(i-1));
bvpfcn = @(eta,y)[y(2);y(3);(4/5)*y(2)^2-(12/5)*y(1)*y(3)-y(4);y(5);-(12/5)*Pr*(y(2)*y(4)+y(1)*y(5))];
if i==1
solinit = bvpinit(etamesh,[0; 0; 0.33206; 1; 0]);
else
solinit = bvpinit(sol_water.x,@(eta)interp1(sol_water.x,(sol_water.y).',eta));
end
sol_water = bvp4c(bvpfcn, bcfcn, solinit);
end
%bvpfcn = @(eta,y)[y(2);y(3);(4/5)*y(2)^2-(12/5)*y(1)*y(3)-y(4);y(5);-(12/5)*Pr_water*(y(2)*y(4)+y(1)*y(5))];
%sol_water = bvp4c(bvpfcn, bcfcn, solinit);
figure(1)
plot(sol_air.x, sol_air.y(2,:), '-o')
hold on
plot(sol_water.x, sol_water.y(2,:), '-o')
hold off
grid on
figure(2)
plot(sol_air.x, sol_air.y(4,:), '-o')
hold on
plot(sol_water.x, sol_water.y(4,:), '-o')
hold off
grid on
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Ordinary Differential Equations dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!