coupled differntial equation using ode45

I'm getting an error while solving this coupled differential equation usually the error is showing issues with vertical concatenation. here's the equation i'm tring to solve with mu0 = exp(-T0) and Boundary conditions as : U(y = -1) = 0 and U(y= 1) = 0 and T0 = 0 at y = -1 and T0 = 1 at y = 1.

here's my code:

% Main script to solve the velocity and temperature profile
clear;
clc;

% Define constants
G = 1;     % Source term (example value)
Na = 1;    % Nusselt number (example value)

% Define the domain for y
y_span = [-1 1];

%% Step 1: Solve Velocity Equation to get U and dU/dy
% The velocity equation is:
% d/dy (mu0 * dU/dy) = G
% Define the velocity equation as a system of two first-order ODEs
function dU = velocity_ode(y, U, mu0, G)
dU = [U(2); (1/mu0) * G];  % U(1) = U, U(2) = dU/dy
end

% Initial conditions for velocity at y = -1
U_initial = [0; 0];  % U = 0 and dU/dy = 0 at y = -1 (you can adjust this)

% Solve the velocity equation using ode45
[y_vel, U_sol] = ode45(@(y, U) velocity_ode(y, U, mu0, G), y_span, U_initial);

Unrecognized function or variable 'mu0'.

Error in solution>@(y,U)velocity_ode(y,U,mu0,G) (line 24)
[y_vel, U_sol] = ode45(@(y, U) velocity_ode(y, U, mu0, G), y_span, U_initial);

Error in odearguments (line 93)
f0 = ode(t0,y0,args{:});

Error in ode45 (line 104)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

% Extract dU/dy from the solution
dU_dy = U_sol(:, 2);  % This is the derivative of U with respect to y

% Plot the velocity profile and its derivative
figure;
subplot(2,1,1);
plot(y_vel, U_sol(:, 1), 'b-', 'LineWidth', 2);  % U(y)
xlabel('y');
ylabel('U(y)');
title('Velocity Profile');

subplot(2,1,2);
plot(y_vel, dU_dy, 'r-', 'LineWidth', 2);  % dU/dy(y)
xlabel('y');
ylabel('dU/dy');
title('Velocity Gradient Profile');

%% Step 2: Solve Temperature Equation using dU/dy from Step 1
% Temperature equation: d^2T0/dy^2 + Na * mu0 * (dU/dy)^2 = 0
% Define the temperature equation as a system of two first-order ODEs
function dT = temperature_ode(y, T, dU_dy, mu0, Na)
dT = zeros(2,1);  % Initialize the output vector

% Interpolate dU/dy from the previously computed solution
dUdy_squared = interp1(y_vel, dU_dy.^2, y, 'linear', 'extrap');

% First equation: dT0/dy = T(2)
dT(1) = T(2);

% Second equation: d^2T0/dy^2 = -Na * mu0 * (dU/dy)^2
dT(2) = -Na * mu0 * dUdy_squared;
end

% Initial conditions for temperature at y = -1
T0_initial = [0; 0];  % T0 = 0 and dT0/dy = 0 at y = -1 (adjust second value if needed)
mu0 = ex(-T0);   % Viscosity (example value)
% Solve the temperature equation using ode45
[y_temp, T_sol] = ode45(@(y, T) temperature_ode(y, T, dU_dy, mu0, Na), y_span, T0_initial);

% Plot the temperature profile
figure;
plot(y_temp, T_sol(:, 1), 'r-', 'LineWidth', 2);  % T0
xlabel('y');
ylabel('T_0(y)');
title('Temperature Profile');

Please help me here, thanks in advance

Answer 1

[Torsten]

G = 1;
Na = 1;

xstart = -1;
xend = 1;
nx = 51;
x = linspace(xstart,xend,nx);
solinit = bvpinit(x, [0;0;1;0]);
sol = bvp4c(@(y,z)bvpfcn(y,z,G,Na), @bcfcn, solinit);
figure(1)
plot(sol.y(1,:), sol.x)
xlabel ('U')
ylabel ('y')

figure(2)
plot(sol.y(3,:), sol.x)
xlabel ('T0')
ylabel ('y')

function dzdy = bvpfcn(y,z,G,Na)
  U     = z(1);
  dUdy  = z(2);
  T0    = z(3);
  dT0dy = z(4);
  dzdy    = zeros(4,1);
  dzdy(1) = dUdy;
  dzdy(2) = dUdy*dT0dy + G*exp(T0);
  dzdy(3) = dT0dy;
  dzdy(4) = -Na*exp(-T0)*dUdy^2;
end

function res = bcfcn(za,zb)
  res = [za(1);zb(1);za(3);zb(3)-1.0];
end

Answer 2

[Shashi Kiran]

Hi @Kushagra Saurabh,

I understand that you are encountering an error while trying to solve the coupled differential equation.

After reviewing your code, here are some suggestions to help resolve the issue.

• Initial mu0 calculation: Update the calculation of mu0 to correctly reflect the initial conditions by using mu0 = exp(-T0_initial(1));. This ensures the viscosity is calculated based on the initial temperature value.

mu0 = exp(-T0_initial(1));   % Viscosity (example value)

• Passing y_vel to temperature_ode: The variable y_vel is used within temperature_ode but isn't passed as an argument. Ensure it's included in the function signature and passed in the ode45 call.

% Solve the temperature equation using ode45
[y_temp, T_sol] = ode45(@(y, T) temperature_ode(y, T, y_vel, dU_dy, mu0, Na), y_span, T0_initial);

• Function Definitions: In MATLAB, functions should be defined at the end of the script or in separate files. Ensure velocity_ode and temperature_ode are properly placed to avoid errors.

Here is the fully executable code incorporating these changes.

% Main script to solve the velocity and temperature profile
clear;
clc;

% Define constants
G = 1;     % Source term (example value)
Na = 1;    % Nusselt number (example value)
mu0 = exp(0);  % Assuming T0 = 0 at y = -1 for initial mu0

% Define the domain for y
y_span = [-1 1];

%% Step 1: Solve Velocity Equation to get U and dU/dy

% Initial conditions for velocity at y = -1
U_initial = [0; 0];  % U = 0 and dU/dy = 0 at y = -1 (you can adjust this)

% Solve the velocity equation using ode45
[y_vel, U_sol] = ode45(@(y, U) velocity_ode(y, U, mu0, G), y_span, U_initial);

% Extract dU/dy from the solution
dU_dy = U_sol(:, 2);  % This is the derivative of U with respect to y

% Plot the velocity profile and its derivative
figure;
subplot(2,1,1);
plot(y_vel, U_sol(:, 1), 'b-', 'LineWidth', 2);  % U(y)
xlabel('y');
ylabel('U(y)');
title('Velocity Profile');

subplot(2,1,2);
plot(y_vel, dU_dy, 'r-', 'LineWidth', 2);  % dU/dy(y)
xlabel('y');
ylabel('dU/dy');
title('Velocity Gradient Profile');

%% Step 2: Solve Temperature Equation using dU/dy from Step 1

% Initial conditions for temperature at y = -1
T0_initial = [0; 0];  % T0 = 0 and dT0/dy = 0 at y = -1 (adjust second value if needed)
mu0 = exp(-T0_initial(1));   % Viscosity (example value)

% Solve the temperature equation using ode45
[y_temp, T_sol] = ode45(@(y, T) temperature_ode(y, T, y_vel, dU_dy, mu0, Na), y_span, T0_initial);

% Plot the temperature profile
figure;
plot(y_temp, T_sol(:, 1), 'r-', 'LineWidth', 2);  % T0
xlabel('y');
ylabel('T_0(y)');
title('Temperature Profile');

%% Functions
% Define the velocity equation as a system of two first-order ODEs
function dU = velocity_ode(y, U, mu0, G)
    % dU = [U(2); (1/mu0) * G];  % U(1) = U, U(2) = dU/dy
    dU = [U(2); (1/mu0) * G];  % U(1) = U, U(2) = dU/dy
end

% Define the temperature equation as a system of two first-order ODEs
function dT = temperature_ode(y, T, y_vel, dU_dy, mu0, Na)
    dT = zeros(2,1);  % Initialize the output vector

    % Interpolate dU/dy from the previously computed solution
    dUdy_squared = interp1(y_vel, dU_dy.^2, y, 'linear', 'extrap');

    % First equation: dT0/dy = T(2)
    dT(1) = T(2);

    % Second equation: d^2T0/dy^2 = -Na * mu0 * (dU/dy)^2
    dT(2) = -Na * mu0 * dUdy_squared;
end

Hope this helps.

[Kushagra Saurabh]

Thanks Shashi for your help but it doesn't gives the solution as shown in this plot

Answer 3

[Torsten]

syms y mu0 G Na U(y) T0(y)
eqn1 = diff(mu0*diff(U,y)) == G;
eqn2 = diff(T0,y,2) + Na*mu0*(diff(U,y))^2 == 0;
conds1 = [U(-1)==0,U(1)==0];
conds2 = [T0(-1)==0,T0(1)==1];
sol = dsolve([eqn1,eqn2],[conds1,conds2])

sol = struct with fields:
    T0: y/2 + (Na*G^2 + 6*mu0)/(12*mu0) - (G^2*Na*y^4)/(12*mu0)
     U: (G*y^2)/(2*mu0) - G/(2*mu0)

If you were told to solve your equations numerically, solve them all together using "bvp4c" and not in two subsequent steps using "ode45". This way, you avoid interpolation of U in the equation for T0 and thus inaccuracies in the solution for T0.

[Kushagra Saurabh]

Thanks Torsten, but I was asked to do using ode45.

[Torsten]

It can't be solved using "ode45" because you have boundary conditions on both endpoints. "ode45" can only be used if both boundary conditions were either given at y = -1 or at y = 1. You have a boundary value problem, not an initial value problem - that's why I suggested to use "bvp4c".

[Kushagra Saurabh]

Ok, so does bvp4c takes into account RK4?

Thanks.

[Torsten]

Ok, so does bvp4c takes into account RK4?

No. Here are the references to the numerical methods used in the code:

References

[1] Shampine, L.F., and J. Kierzenka. "A BVP Solver based on residual control and the MATLAB PSE." ACM Trans. Math. Softw. Vol. 27, Number 3, 2001, pp. 299–316.

[2] Shampine, L.F., M.W. Reichelt, and J. Kierzenka. "Solving Boundary Value Problems for Ordinary Differential Equations in MATLAB with bvp4c." MATLAB File Exchange, 2004.

[Kushagra Saurabh]

Thanks Torsten. l'll have a look.