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

Question

sajjad 2024-9-21

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2154545-how-to-solve-the-system-of-time-dependent-coupled-pde-s

评论： sajjad 2024-10-9

在 MATLAB Online 中打开

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

% 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

% 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 个评论
显示 2更早的评论隐藏 2更早的评论

sajjad 2024-9-22

移动：Torsten 2024-9-22

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

sajjad 2024-9-22

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

请先登录，再进行评论。

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

Answer 1

Torsten 2024-9-22

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2154545-how-to-solve-the-system-of-time-dependent-coupled-pde-s#answer_1520800

编辑：Torsten 2024-9-23

在 MATLAB Online 中打开

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

23 个评论
显示 21更早的评论隐藏 21更早的评论

Torsten 2024-9-23

编辑：Torsten 2024-9-23

在 MATLAB Online 中打开

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

Torsten 2024-9-24

编辑：Torsten 2024-9-24

在 MATLAB Online 中打开

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 2024-9-24

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.

Torsten 2024-9-26

编辑：Torsten 2024-9-26

在 MATLAB Online 中打开

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 2024-9-27

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

Torsten 2024-9-27

编辑：Torsten 2024-9-27

在 MATLAB Online 中打开

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

Torsten 2024-10-8

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 2024-10-9

yes right equations should be known.

请先登录，再进行评论。

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

4 个评论
显示 2更早的评论隐藏 2更早的评论

回答（1 个）

23 个评论
显示 21更早的评论隐藏 21更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

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

4 个评论 显示 2更早的评论隐藏 2更早的评论

回答（1 个）

23 个评论 显示 21更早的评论隐藏 21更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

4 个评论
显示 2更早的评论隐藏 2更早的评论

23 个评论
显示 21更早的评论隐藏 21更早的评论