Solving the function ODE45
1 view (last 30 days)
Show older comments
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 Comment
Torsten
on 13 Jun 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).
Answers (1)
Nipun
on 13 Jun 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 Comment
John D'Errico
on 13 Jun 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.
See Also
Categories
Find more on Eigenvalue Problems in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!