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

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

[Torsten]

I can't understand the update scheme in time that you use for F and G. Where do you solve for fvar and gvar as solutions for time (i+1)*dt ?

Why don't you follow the numerical procedure (3.1) - (3.8) explained in the article ?

[sajjad]

Actually i need just the output of F and G so i used FDM. The method in the article i was not understanding. I possible help to solve the system with the method in the article. Thank you.

[sajjad]

First i have to draw fig 10a and then Fig 8. Then I will use the outputs for my further research work.

[sajjad]

In equation 3.9 he is describing the boundary conditions for G as well.

Answer 1

[Torsten]

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:

http://www.unige.ch/~hairer/software.html

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

[sajjad]

well. If possible put the whole code in one window.? I tried to run it but invain errors are occuring. In this code can be Renolds and Delta values be changed? I think it shuold be general code, that can be changed for every fig and differnet parametrs.

[Torsten]

If you get an error message, include it here in its full length. Which MATLAB release do you use ?

delta and Reynolds can be changed in the script part.

[sajjad]

It is runing online but not on system. >> mainScript

Error: File: mainScript.m Line: 28 Column: 1

Function definitions are not permitted in this context.

[Torsten]

Save the function in a separate file and name it fun.m. Delete the function from the main script and try again.

With "function" I mean this part:

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

[sajjad]

well.

Now its running sucessfully but slowly, but a little differnce in article the range of y-axis in fig10a is almot 0 to 0.2, but in code it is little bit diffrent.

I am confuse that our system is of PDE's how we will use ODE23t? Can our system, intial and boundary conditions in code be explaind? This code is for fig10a, i.e the output of G(1/4,t) y=1/4 and t from 8100 to 8200, if want to see the outputs of G(1/4,t) how i can? One thing more if want find the outputs of G(y,t) at spscific value of y i.e in Fig 2 (a,b) maxima and minima of α(t)=G(1,t) at y=1 and their corresponding return maps at delta =0.3 and Relnolds=900, how it will plot.?

[sajjad]

This is equation number (2.4 b).

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);

i m not understaning this. In (2.4b) there is no derivative w.r.t t.

[Torsten]

I am confuse that our system is of PDE's how we will use ODE23t? Can our system, intial and boundary conditions in code be explaind?

The code applies the method-of-lines to your PDE system.

After discretizing your equations in space, you get a mixture of ordinary differential equations and algebraic equations in the grid points. For F, you get algebraic equations because there is no time derivative in the elliptic equation for F, for G (except for the first and for the last grid point) you get differential equations in the grid points because of the parabolic equation for G.

ode15s and ode23t solve systems of differential-algebraic equations of the form

M*y' = f(t,x)

where M is the so called mass matrix. The zero rows of M correspond to the algebraic equations, the rows with a 1 correspond to the differential equations.

If you look at how M is defined in the code, you will see that M is arranged such that the equations solved become

F(1) = 0;
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) = 0
F(nx) + Hdot(t) = 0
G(1) = 0;
dG/dt(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);
G(nx) - (d2Fdx2(nx)+dFdx(nx)/x(nx)-F(nx)/x(nx)^2)/(-H(t)^2) = 0

I admit that the name "dFdt" is misleading here because dFdt are not time-derivatives, but algebraic expressions that should be solved to dFdt == 0. Same for dGdt(1) and dGdt(nx).

This code is for fig10a, i.e the output of G(1/4,t) y=1/4 and t from 8100 to 8200, if want to see the outputs of G(1/4,t) how i can?

Replace

idx = T >= 8100 

by

idx = T >= 0

One thing more if want find the outputs of G(y,t) at spscific value of y i.e in Fig 2 (a,b) maxima and minima of α(t)=G(1,t) at y=1 and their corresponding return maps at delta =0.3 and Relnolds=900, how it will plot.?

[pks,locs] = findpeaks(G(:,nx))
plot(T(locs),pks,'o')

plots the maxima of G(eta=1) over time

[pks,locs] = findpeaks(-G(:,nx))
pks = -pks
plot(T(locs),pks,'o')

plots the minima of G(eta=1) over time.

If you don't manage the plotting part and you want to further work on such a complicated problem as the one from above, I suggest you first learn the basics of MATLAB and pass MATLAB Onramp, an introductory course in MATLAB free of costs:

https://matlabacademy.mathworks.com/details/matlab-onramp/gettingstarted

A general remark:

The presets for ode23t are a relative tolerance of 1e-3 and an absolute tolerance of 1e-6. For such a strongly dynamic system as yours it seems advisable to strengthen these values in the odeset().

Also playing with the number of discretization points nx is a necessary step to see how the solutions behave.

[sajjad]

Perfect, i understood that points. Thank you.

Last point i have sucessfully plotted the maxima and minima of G(1,t) but polt is not simialr to fig2a and fig2b, for this i am sending the modify code.

Please guide me how to plot return maps of maxima and mimina?

Thank you.

clear all;

close all;

clc;

nx = 101;

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 = 900;

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 >= 0;

figure(1)

plot(T(idx),G(:,nx));

[pks,locs] = findpeaks(G(:,nx));

figure(2)

plot(T(locs),pks,'o')

[pks,locs] = findpeaks(-G(:,nx))

pks = -pks

figure(3)

plot(T(locs),pks,'o')

Sir, and also print the data (inputs on which the maps are made) of return maps.

[sajjad]

[Torsten]

I got the same plot with this code and with a different integrator as well. So I suspect that the plot in the article is wrong. Maybe the scale on the time axis should be in units of 100 instead of units 1000.

Your included plot of the minima does not correspond to what comes from the code: The values on the y-axis should be from -58 to -54, not from 54 to 58.

Sir, and also print the data (inputs on which the maps are made) of return maps.

Sorry, no. Making plots is boring.

[Torsten]

I found out that the boundary between chaotic and regular behaviour strongly depends on the error tolerances for the integrator and the mesh size (value for nx) in eta-direction. Having this in mind, I wonder if numerical computations with limited accuracy of quasi-chaotic systems make sense at all.

[sajjad]

Okay.

[sajjad]

(1)

The most important thing i want to discuss, the relation between number of iteratations , time range and dx = (xend-xstart)/(nx-1), here we can see that dx=0.01, for this if i need the out put for t from 0 to 8000 then what shoulb be the the number of iterations iteractions? i.e in fig 2a. tis from 0 to 8000. so for this t_end will be?

(2)

In the expression dx = (xend-xstart)/(nx-1), 1 will be always fix?

(3)

In the expresion plot(T(idx),G(:,nx)); this is for G(1,t), and if i ned to plot the output of G(1/2,t) for any range of t, then what will be chnage in plot(T(idx),G(:,nx))?

(4)

How will nx can be change?

[Torsten]

nx is the number of grid points of the mesh between eta = 0 and eta = 1 (eta is called x in the code).

tend is the time up to which the solver shall integrate.

dx is computed from nx and xstart and xend and is the length of a single cell in the eta-mesh. Given xstart, xend and nx, you can easily deduce that it is always given by dx = (xend-xstart)/(nx-1).

In the expresion plot(T(idx),G(:,nx)); this is for G(1,t), and if i ned to plot the output of G(1/2,t) for any range of t, then what will be chnage in plot(T(idx),G(:,nx))?

If you have nx grid points equally spaced between 0 and 1, can you compute which grid point represents eta = 1/2 ? It's grid point number (nx+1)/2. So you should arrange nx such that nx+1 is divisible by 2.

For eta = 1/4, it's grid point number (nx+3)/4, I guess. So in this case, you should arrange nx such that nx+3 is divisible by 4 (as 26 in the example code with nx = 101).

And the command must be

plot(T(idx),G(idx,nx))

, not

plot(T(idx),G(:,nx))

How will nx can be change?

It's defined in the first line of the code.

[sajjad]

Here in code figure 2 a, b (maxima and minima of G(1,t) are not scatrering) after 1000 time. I need the results t from 0 to 8000 figs has been attached. Thank you.

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-2);

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 = 900;

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 >= 0;

figure(1)

plot(T(idx),G(:,nx));

xlabel('Time (t)');

ylabel('G(1,t)');

%title('G(1,t)');

grid on;

% Assuming G and T have already been defined and filtered for t <= 800

G1 = G (: , 101); % 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', 4); % 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', 4); % 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

[sajjad]

How to drwa/contruct the geometry of the problem of the article?

[Torsten]

If you don't get the results for F and G as in the article, you can't get the plots for maxima, minima and return maps as in the article.

In my opinion, this is the equivalent to what is plotted in fig.2.

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.3;
Reynolds = 900;
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);
[pks_max,locs_max] = findpeaks(G(:,nx));
[pks_min,locs_min] = findpeaks(-G(:,nx));
pks_min = -pks_min;
figure(1)
plot(T(locs_max),pks_max,'.')

figure(2)
plot(pks_max(1:end-1),pks_max(2:end),'.')

figure(3)
plot(T(locs_min),pks_min,'.')

figure(4)
plot(pks_min(1:end-1),pks_min(2:end),'.')

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

[sajjad]

Sir i am wondering that from this code we are getting the same figure 10 (a) in the article, but why other results are not matching?

[Torsten]

We don't have plots as in fig. 2a and b for the conditions delta = 0.2 and R = 2050.

The time span of the chaotic regime for the conditions in fig. 2a and b is from 0 to 400 appr. with the code from above whereas it is from 0 to 4000 appr. in the author's publication. So I'd say that the discrepancy is large in all respects.

As I wrote several times now, you must vary the tolerances and the number of mesh points. The aim is to get stable solutions that don't change if tolerances are chosen even more stringent and if the number of mesh points is further increased. And as I also noted, according to my simulations, this does not happen. My conclusion would be that numerical simulations for this system in the moderate to high Reynolds number regime are useless.

[sajjad]

Thank you. I have almost draw the figures same like article.

[sajjad]

The problem/system discussed earliar was the unsteady motion of a viscous incompressible fluid inside a cylindrical tube whose radius is changing in a prescribed manner.

Can we modify our code or system for Non-Newtonian fluid case or some other type of flow?

[Torsten]

Can we modify our code or system for Non-Newtonian fluid case or some other type of flow?

Maybe - if you know the new equations :-)

[sajjad]

yes right equations should be known.