Solving the function ODE45

1 vue (au cours des 30 derniers jours)
JungYeon Han
JungYeon Han le 13 Juin 2024
Commenté : John D'Errico le 13 Juin 2024
I have to solve the following equation with the initial condition of C_{11}=1 and C_{12}=0.
Large D means material derivative.
Following is the code, and I am not sure how I could run the code.
function solve_pde()
% Parameters
Gamma = 4.0; % Circulation
nu = 0.01; % Kinematic viscosity
x = linspace(-2, 2, 100); % x-coordinates
y = linspace(-2, 2, 100); % y-coordinates
[X, Y] = meshgrid(x, y);
tspan = [0, 1]; % Time span for the simulation
% Initial conditions for the conformation tensor components
C_11_init = ones(size(X));
C_12_init = zeros(size(X));
% Combine initial conditions into a single vector for the ODE solver
C0 = [C_11_init(:); C_12_init(:)];
% Solve the system of ODEs
[t, C] = ode45(@(t, C) conformation_tensor_pde(t, C, X, Y, Gamma, nu), tspan, C0);
% Extract the results for each component of the conformation tensor
num_points = numel(X);
C_11 = reshape(C(end, 1:num_points), size(X));
C_12 = reshape(C(end, num_points+1:end), size(X));
% Plot the results
figure
surf(X, Y, C_11)
xlabel('X');
ylabel('Y');
zlabel('C_{11}');
title('Evolution of C_{11}');
end
function dCdt = conformation_tensor_pde(t, C, X, Y, Gamma, nu)
% Reshape C into the conformation tensor components
num_points = numel(X);
C_11 = reshape(C(1:num_points), size(X));
C_12 = reshape(C(num_points+1:end), size(X));
% Velocity field components in Cartesian coordinates
r = sqrt(X.^2 + Y.^2);
u = -Gamma * Y ./ (2 * pi * (X.^2 + Y.^2)) .* (1 - exp(-r.^2 / (4 * nu * t)));
v = Gamma * X ./ (2 * pi * (X.^2 + Y.^2)) .* (1 - exp(-r.^2 / (4 * nu * t)));
% Velocity gradients in Cartesian coordinates
dudx = Gamma * Y .* (2 * X .* (1 - exp(-r.^2 / (4 * nu * t))) + X .* exp(-r.^2 / (4 * nu * t)) / (2 * nu * t)) ./ (2 * pi * (X.^2 + Y.^2).^2);
dudy = -Gamma * (1 - exp(-r.^2 / (4 * nu * t))) ./ (2 * pi * (X.^2 + Y.^2)) + Gamma * Y.^2 .* (2 * (1 - exp(-r.^2 / (4 * nu * t))) - exp(-r.^2 / (4 * nu * t)) / (2 * nu * t)) ./ (2 * pi * (X.^2 + Y.^2).^2);
% Finite differences for spatial derivatives
[dC11dx, dC11dy] = gradient(C_11, X(1, :), Y(:, 1));
% Compute the material derivative of the conformation tensor component C_11
dC11dt = 2 * (C_11 .* dudx + C_12 .* dudy) - u .* dC11dx - v .* dC11dy;
% For completeness, compute the derivative of C_12 (if necessary for other equations)
dC12dt = zeros(size(C_12)); % Adjust based on the specific problem requirements
% Pack the derivatives into a column vector
dCdt = [dC11dt(:); dC12dt(:)];
end
  1 commentaire
Torsten
Torsten le 13 Juin 2024
Are u and v given function ?
You only have one partial differential equation for C_11. Where is the equation for C_12 ?
Further, assuming u and v are positive, you need boundary conditions at the lines (xmin,y) and (x,ymin).

Connectez-vous pour commenter.

Réponses (1)

Nipun
Nipun le 13 Juin 2024
Hi JungYeon,
I understand that you intend to solve the given ordinary differential equation using MATLAB and require help with running the given code.
Based on the given code and upon close inspection, I can verify that the given code corresponds to the solution for the given ODE.
In order to call your solution, type the following in the MATLAB command window:
>> solve_pde();
This command will call the function "solve_pde", which in turn uses "conformation_tensor_pde" to solve the equation. Additionally, a plot will be generated at the end of the function.
For more information on calling functions in MATLAB, refer to the following MathWorks documentation: https://www.mathworks.com/help/matlab/ref/function.html
Hope this helps.
Regards,
Nipun
  1 commentaire
John D'Errico
John D'Errico le 13 Juin 2024
That is NOT an ordernary differential equation. It is a PARTIAL differential equation. There are derivatives with respect to three variables, t, x, and y.

Connectez-vous pour commenter.

Tags

Community Treasure Hunt

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

Start Hunting!

Translated by