Unable to meet integration tolerances
3 views (last 30 days)
Show older comments
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 Comments
Accepted Answer
Torsten
on 4 Feb 2024
Edited: Torsten
on 4 Feb 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 Comments
Torsten
on 4 Feb 2024
Edited: Torsten
on 4 Feb 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
More Answers (0)
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!