coupled differntial equation using ode45

9 次查看（过去 30 天）
Kushagra Saurabh
Kushagra Saurabh 2024-9-27
评论： Kushagra Saurabh 2024-9-30
Ran in:
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

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

采纳的回答

Torsten
Torsten 2024-9-27
编辑：Torsten 2024-9-27
Ran in:
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

更多回答（2 个）

Shashi Kiran
Shashi Kiran 2024-9-27
Ran in:
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 个评论
Kushagra Saurabh
Kushagra Saurabh 2024-9-27
Thanks Shashi for your help but it doesn't gives the solution as shown in this plot

请先登录，再进行评论。

Torsten
Torsten 2024-9-27
编辑：Torsten 2024-9-27
Ran in:
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更早的评论
Torsten
Torsten 2024-9-28
编辑：Torsten 2024-9-28
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by