Filter löschen
Filter löschen

Unable to meet integration tolerances

3 Ansichten (letzte 30 Tage)
Panda05
Panda05 am 3 Feb. 2024
Bearbeitet: Torsten am 4 Feb. 2024
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
Torsten
Torsten am 3 Feb. 2024
Bearbeitet: Torsten am 3 Feb. 2024
Why do you use ode45 and not bvp4c/bvp5c, the solvers for boundary value problems in MATLAB ?
And what makes you think that 0 < f''(0) < 1 ?
Panda05
Panda05 am 3 Feb. 2024
I have never used bvp4c before unfortunately . when i tried to think of a way to solve this problem , i used 4th order Runge Kuntta as its the one I was familiar with . But if this can help solve the problem , I'd be greatful to show me how

Melden Sie sich an, um zu kommentieren.

Akzeptierte Antwort

Torsten
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
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
Panda05
Panda05 am 4 Feb. 2024
Thank you so much . This is much better ! You put an end to the powerful rage that was building up inside me ;).

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (0)

Tags

Community Treasure Hunt

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

Start Hunting!

Translated by