Unable to meet integration tolerances
3 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
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 Kommentare
Akzeptierte Antwort
Torsten
am 4 Feb. 2024
Bearbeitet: Torsten
am 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 Kommentare
Torsten
am 4 Feb. 2024
Bearbeitet: Torsten
am 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
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Ordinary Differential Equations finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!