Unable to meet integration tolerances

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

[Torsten]

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]

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

Answer 1

[Torsten]

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

[Panda05]

Hello , thank you so much for the solution . The thing is the Pr value should not be changed since its a constant in the original question I'm solving .

I tried to use python, it is working , however , i wanted to make an equivalent code here in Matlab and couldnt come up with a good code

import numpy as np

from scipy.integrate import solve_ivp

# 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

# Define the system of ODEs

def odes_system(eta, y, Pr):

f, X, Y, theta, Z = y

dfdeta = X

dXdeta = Y

dYdeta = (4/5)*X**2 - (12/5)*f*Y - theta # f''' = ...

dthetadeta = Z

dZdeta = -(12/5)*Pr*(X*theta + f*Z)

return [dfdeta, dXdeta, dYdeta, dthetadeta, dZdeta]

# 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

# Function to execute the Runge-Kutta integration and adjust f''(0) iteratively

def solve_with_adjustment(Pr, y0, eta_end, delta_eta):

# 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

eta_span = [0, eta_end]

for i in range(max_iterations):

# Solve the ODEs

y0[2] = (f_double_prime_0_lower + f_double_prime_0_upper) / 2 # Update the guess for f''(0)

sol = solve_ivp(odes_system, eta_span, y0, args=(Pr,), method='RK45', max_step=delta_eta)

# Check the boundary condition at infinity

if abs(sol.y[1, -1]) < tolerance: # If f'(eta_end) is close enough to 0

return sol

elif sol.y[1, -1] > 0: # f'(eta_end) is too high, decrease f''(0)

f_double_prime_0_upper = y0[2]

else: # f'(eta_end) is too low, increase f''(0)

f_double_prime_0_lower = y0[2]

raise ValueError("Solution did not converge within the maximum number of iterations")

# 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

plt.figure(figsize=(12, 6))

# Velocity profile

plt.subplot(1, 2, 1)

plt.plot(sol_air.t, sol_air.y[1, :], label=f'Air (Pr={Pr_air})')

plt.plot(sol_water.t, sol_water.y[1, :], label=f'Water (Pr={Pr_water})')

plt.title('Velocity $f\'(\eta)$')

plt.xlabel('$\eta$')

plt.ylabel("$f'(\eta)$")

plt.legend()

# Temperature profile

plt.subplot(1, 2, 2)

plt.plot(sol_air.t, sol_air.y[3, :], label=f'Air (Pr={Pr_air})')

plt.plot(sol_water.t, sol_water.y[3, :], label=f'Water (Pr={Pr_water})')

plt.title(' Temperature $\\theta(\eta)$')

plt.xlabel('$\eta$')

plt.ylabel('$\\theta(\eta)$')

plt.legend()

plt.tight_layout()

plt.show()

[Torsten]

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:

https://uk.mathworks.com/matlabcentral/fileexchange/64968-rk45-dydt-trange-yinit-tol

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]

Thank you so much . This is much better ! You put an end to the powerful rage that was building up inside me ;).