coupled differntial equation using ode45
4 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Kushagra Saurabh
el 27 de Sept. de 2024
Comentada: Kushagra Saurabh
el 30 de Sept. de 2024
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);
% 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 comentarios
Respuesta aceptada
Torsten
el 27 de Sept. de 2024
Editada: Torsten
el 27 de Sept. de 2024
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 comentarios
Más respuestas (2)
Shashi Kiran
el 27 de Sept. de 2024
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.
Torsten
el 27 de Sept. de 2024
Editada: Torsten
el 27 de Sept. de 2024
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])
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 comentarios
Torsten
el 28 de Sept. de 2024
Editada: Torsten
el 28 de Sept. de 2024
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.
Ver también
Categorías
Más información sobre Ordinary Differential Equations en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!