How to solve the system of time dependent coupled PDE's?

4 vues (au cours des 30 derniers jours)
sajjad
sajjad le 21 Sep 2024
The system this paper (DOI: 10.1017/S0022112003003835) Thank you.
Fig10a(0.2,2050)
Iteration: 81000 | G(1/4, t): NaN Iteration: 82000 | G(1/4, t): NaN
function Fig10a(delta, Reynolds)
% Main function to solve for F and G and plot G(1/4, t)
% Initialize variables
R = Reynolds;
dt = 1/10;
nmax = 82001;
dy = 1/100;
yTarget = 1/4;
gValues = [];
Delta = delta;
H = @(t) 1 + Delta * cos(2 * t);
dH = @(t) -2 * Delta * sin(2 * t);
% Create the grid and differentiation matrices
ygrid = 0:dy:1;
ny = length(ygrid);
% Finite difference differentiation matrices (for dy1 and dy2)
dy2 = fdcoeffFDM2(ny, dy); % Second-order finite difference matrix
dy1 = fdcoeffFDM1(ny, dy); % First-order finite difference matrix
% Initialize variables for F and G
fvar0 = zeros(ny, nmax);
gvar0 = zeros(ny, nmax);
% Time-stepping loop
for i = 1:nmax-1
% Update variables for F and G at the current time step
fvar = fvar0(:, i);
gvar = gvar0(:, i);
% Define the equations
eqf = dy2 * fvar + (dy1 * fvar)./ygrid' - fvar./(ygrid'.^2) + ...
gvar * H((i + 1) * dt)^2;
eqg = (gvar - gvar0(:, i)) / dt - ...
dH((i + 1) * dt) / H((i + 1) * dt) * (dy1 * gvar)./ygrid' - ...
(dy1 * gvar) .* fvar / H((i + 1) * dt) + ...
1/H((i + 1) * dt) * (dy1 * fvar) .* gvar + ...
2./(ygrid' * H((i + 1) * dt)) .* fvar .* gvar - ...
(1/(R * H((i + 1) * dt)^2)) * (dy2 * gvar + dy1 * gvar./ygrid' - gvar./(ygrid'.^2));
% Apply boundary conditions
eqf(1) = fvar(1); % F(0, t) = 0 at y = 0
eqf(end) = fvar(end) + dH((i + 1) * dt); % F(1, t) = -H'(t) at y = 1
eqg(1) = gvar(1); % G(0, t) = 0 at t = 0
eqg(end) = dy1(end, :) * fvar - dH((i + 1) * dt); % F'(1, t) = H'(t) at y = 1
% Combine eqf and eqg into a single system
eqns = [eqf; eqg];
% Solve the system using backslash operator
sol = eqns; % Since we are already evaluating eqns as the result
% Check that the solution vector matches the expected size
if length(sol) ~= 2 * ny
error('The number of equations does not match the number of unknowns');
end
% Update variables for the next step
fvar0(:, i+1) = sol(1:ny);
gvar0(:, i+1) = sol(ny+1:end);
% Interpolate G(y, t) and store G(1/4, t)
if i >= 81000 && i <= 82001
gInterp = interp1(ygrid, gvar0(:, i+1), yTarget);
gValues = [gValues; i, gInterp];
end
% Display debugging information every 1000 iterations
if mod(i, 1000) == 0 & i>=81000
disp(['Iteration: ', num2str(i), ' | G(1/4, t): ', num2str(gInterp)]);
end
end
% Plot the values of G(1/4, t)
figure;
if isempty(gValues)
disp('No values of G(1/4, t) were recorded.');
else
plot(gValues(:,1), gValues(:,2), 'r-', 'LineWidth', 2);
grid on;
xlabel('t');
ylabel('G(1/4, t)');
title('G(1/4, t) from t = 8100 to t = 8200');
end
end
% Auxiliary function for second-order finite difference matrix
function D2 = fdcoeffFDM2(ny, dy)
e = ones(ny, 1);
D2 = spdiags([e -2*e e], -1:1, ny, ny) / (dy^2);
D2(1, :) = 0; D2(end, :) = 0; % Apply boundary conditions
end
% Auxiliary function for first-order finite difference matrix
function D1 = fdcoeffFDM1(ny, dy)
e = ones(ny, 1);
D1 = spdiags([-e e], [-1 1], ny, ny) / (2*dy);
D1(1, :) = 0; D1(end, :) = 0; % Apply boundary conditions
end
  4 commentaires
sajjad
sajjad le 22 Sep 2024
Déplacé(e) : Torsten le 22 Sep 2024
First i have to draw fig 10a and then Fig 8. Then I will use the outputs for my further research work.
sajjad
sajjad le 22 Sep 2024
In equation 3.9 he is describing the boundary conditions for G as well.

Connectez-vous pour commenter.

Réponses (1)

Torsten
Torsten le 22 Sep 2024
Modifié(e) : Torsten le 23 Sep 2024
ode23t seems to work for high Reynolds numbers. In case ode23t fails (e.g. for the low Reynolds number regime), I recommend the solver "radau" for MATLAB which seems to works for the complete Reynolds number regime, but is much slower than ode23t.
It can be downloaded from:
nx = 101;
xstart = 0.0;
xend = 1.0;
x = linspace(xstart,xend,nx).';
dx = (xend-xstart)/(nx-1);
tstart = 0;
tend = 8200;
tspan = [tstart,tend];
G0 = zeros(nx,1);
F0 = zeros(nx,1);
y0 = [G0;F0];
delta = 0.2;
Reynolds = 2050;
H = @(t) 1 + delta*cos(2*t);
Hdot = @(t) -delta*2*sin(2*t);
M = [eye(nx),zeros(nx);zeros(nx,2*nx)];
M(1,1) = 0;
M(nx,nx) = 0;
options = odeset('Mass',M);
[T,Y] = ode23t(@(t,y)fun(t,y,x,nx,dx,delta,Reynolds,H,Hdot),tspan,y0,options);
G = Y(:,1:nx);
F = Y(:,nx+1:2*nx);
idx = T >= 8100;
figure(1)
plot(T(idx),G(idx,26))
figure(2)
plot(T(idx),F(idx,26))
function dydt = fun(t,y,x,nx,dx,delta,Reynolds,H,Hdot)
%persistent iter
%if isempty(iter)
% iter = 0;
%end
%iter = iter + 1;
%if mod(iter,1000)==0
% t
% iter = 0;
%end
G = y(1:nx);
F = y(nx+1:2*nx);
% F
% Compute spatial derivatives
dFdx = zeros(nx,1);
d2Fdx2 = zeros(nx,1);
dFdx(2:nx-1) = (F(3:nx)-F(1:nx-2))/(2*dx);
dFdx(nx) = Hdot(t);
d2Fdx2(2:nx-1) = (F(3:nx)-2*F(2:nx-1)+F(1:nx-2))/dx^2;
%d2Fdx2(nx) = (dFdx(nx)-dFdx(nx-1))/dx;
Fnxp1 = F(nx-1) + 2*dx*dFdx(nx);
d2Fdx2(nx) = (Fnxp1-2*F(nx)+F(nx-1))/dx^2;
% Compute temporal derivatives
dFdt = zeros(nx,1);
dFdt(1) = F(1);
dFdt(2:nx-1) = d2Fdx2(2:nx-1)+dFdx(2:nx-1)./x(2:nx-1)-...
F(2:nx-1)./x(2:nx-1).^2+H(t)^2*G(2:nx-1);
dFdt(nx) = F(nx) + Hdot(t);
% G
% Compute spatial derivatives
dGdx = zeros(nx,1);
d2Gdx2 = zeros(nx,1);
dGdx(2:nx-1) = (G(3:nx)-G(1:nx-2))/(2*dx);
d2Gdx2(2:nx-1) = (G(3:nx)-2*G(2:nx-1)+G(1:nx-2))/dx^2;
% Compute temporal derivatives
dGdt = zeros(nx,1);
dGdt(1) = G(1);
dGdt(2:nx-1) = Hdot(t)/H(t).*x(2:nx-1).*dGdx(2:nx-1)+ ...
1/H(t).*F(2:nx-1).*dGdx(2:nx-1)- ...
1/H(t).*dFdx(2:nx-1).*G(2:nx-1)- ...
2.*F(2:nx-1).*G(2:nx-1)./(H(t)*x(2:nx-1))+ ...
(d2Gdx2(2:nx-1)+dGdx(2:nx-1)./x(2:nx-1)-G(2:nx-1)./x(2:nx-1).^2)./ ...
(H(t)^2*Reynolds);
dGdt(nx) = G(nx) - (d2Fdx2(nx)+dFdx(nx)/x(nx)-F(nx)/x(nx)^2)/(-H(t)^2);
%Taken from the article, should be the same as set
%dGdt(nx) = G(nx) + 2/H(t)^2*((F(nx-1)+Hdot(t)*(1+dx))/dx^2 + Hdot(t));
dydt = [dGdt;dFdt];
end
  25 commentaires
sajjad
sajjad le 16 Déc 2024
Modifié(e) : sajjad le 16 Déc 2024
In the below code i want to increase the number of outputs till 65563 for further work. How is it posible.? This is the code for fig 3 in the previous papper graph of G(eeta,t), at eeta = 1. i am not understating the how to increase or decrea the discretization?
sajjad
sajjad le 16 Déc 2024
Modifié(e) : Walter Roberson le 15 Mar 2025
clear all;
close all;
clc;
nx = 101; %The number of points or grid cells used to discretize the spetial domain
xstart = 0.0;
xend = 1.0;
x = linspace(xstart,xend,nx).';
dx = (xend-xstart)/(nx-1);
tstart = 0;
tend =8000;
tspan = [tstart,tend];
G0 = zeros(nx,1); % corresponding to F(etta,0)=0
F0 = zeros(nx,1); % corresponding to G(etta,0)=0
y0 = [G0;F0];
delta = 0.3;
Reynolds = 1000;
H = @(t) 1 + delta*cos(2*t);
Hdot = @(t) -delta*2*sin(2*t);
M = [eye(nx),zeros(nx);zeros(nx,2*nx)];
M(1,1) = 0;
M(nx,nx) = 0;
options = odeset('Mass',M);
[T,Y] = ode23t(@(t,y)funn(t,y,x,nx,dx,delta,Reynolds,H,Hdot),tspan,y0,options);
G = Y(:,1:nx);
F = Y(:,nx+1:2*nx);
idx = T >= 0;
figure(1)
plot(T(idx),G(:,nx));
% Assuming G and T have already been defined and filtered for t <= 800
G1 = G (: , nx); % Values of G at eta = 1
t_range = T (idx); % Corresponding time points
% Find maxima
[maxima, locs_max] = findpeaks (G1 (idx)); % Find peaks and their locations
max_times = t_range (locs_max); % Times at which maxima occur
% Find minima
[minima, locs_min] = findpeaks (-G1 (idx)); % Negate to find minima
min_times = t_range (locs_min); % Times at which minima occur
% Plot Maxima as points only
figure(2);
plot(max_times, maxima, 'ro', 'MarkerFaceColor', 'r','MarkerSize', 3) % Plot maxima as red points
xlabel('Time (t)');
ylabel('Maxima of G(1,t)');
title('Maxima of G(1,t)');
grid on;
% Plot Minima as points only
figure(3);
plot (min_times, -minima, 'go', 'MarkerFaceColor', 'g','MarkerSize', 3) % Plot minima as green points
xlabel (' Time (t)');
ylabel (' Minima of G (1, t)');
title (' Minima of G (1, t)');
grid on;
[pks1,locs] = findpeaks(G(:,nx));
figure(4)
plot(T(locs),pks1,'o')
[pks2,locs] = findpeaks(-G(:,nx));
%pks = -pks
figure(5)
plot(T(locs),pks2,'o')
% Create return map for maxima: plot max(i+1) vs max(i)
figure (6)
plot(maxima(1:end-1), maxima(2:end), 'ro', 'MarkerFaceColor', 'r','MarkerSize', 3); % Return map for maxima
xlabel('Max(i)');
ylabel('Max(i+1)');
title('Return Map of Maxima');
grid on;
% Display data for return map
return_map_maxima = [maxima(1:end-1), maxima(2:end)]; % Prepare data for return map
disp('Return Map Data for Maxima (Max(i), Max(i+1)):');
disp(return_map_maxima); % Display the consecutive maxima pairs in the Command Window
% Create return map for minima: plot min(i+1) vs min(i)
figure(7)
plot(-minima(1:end-1), -minima(2:end), 'go', 'MarkerFaceColor', 'g','MarkerSize', 3); % Return map for minima
xlabel('Min(i)');
ylabel('Min(i+1)');
title('Return Map of Minima');
grid on;
% Display data for return map
return_map_minima = [minima(1:end-1), minima(2:end)]; % Prepare data for return map
disp('Return Map Data for Maxima (Max(i), Max(i+1)):');
disp(return_map_minima); % Display the consecutive maxima pairs in the Command Window

Connectez-vous pour commenter.

Catégories

En savoir plus sur Numerical Integration and Differential Equations dans Help Center et File Exchange

Produits


Version

R2013a

Community Treasure Hunt

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

Start Hunting!

Translated by