coupled differntial equation using ode45

Question

Kushagra Saurabh 2024-9-27，6:49

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2155800-coupled-differntial-equation-using-ode45

评论： Kushagra Saurabh 2024-9-30，4:53

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2024-9-27，12:42

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2155800-coupled-differntial-equation-using-ode45#answer_1523735

编辑：Torsten 2024-9-27，17:03

在 MATLAB Online 中打开

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Answer 2

Shashi Kiran 2024-9-27，7:46

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2155800-coupled-differntial-equation-using-ode45#answer_1523505

在 MATLAB Online 中打开

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.

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Kushagra Saurabh 2024-9-27，8:33

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

请先登录，再进行评论。

Answer 3

Torsten 2024-9-27，8:34

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2155800-coupled-differntial-equation-using-ode45#answer_1523530

编辑：Torsten 2024-9-27，8:45

在 MATLAB Online 中打开

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.

5 个评论
显示 3更早的评论隐藏 3更早的评论

Torsten 2024-9-28，10:10

编辑：Torsten 2024-9-28，14:51

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 2024-9-30，4:53