Implicit finite difference for 5 simultaneous pde and Newton-Raphson method

6 views (last 30 days)
Hello,
I am trying to solve my system with 5 nonlinear pde with 5 unknown functions using implicit finite difference method. At the same time, the code uses Newton-Raphson iteration for gap1_w+gap2_w=1. I have coded the problem as shown below
%--------------------------------------------------------------------------
Nt = 3; % Total time step number
Nx = 101; % Total grid
dt = 0.0001; % Time step
dx = 1/(Nx-1); % Grid size
xplot = 1:dx:2; % Record the x scale for plots
r = dt/dx; % define the ratio r
iplot = 1; % Counter used to count plots
tplot(1) = 0; % Initial time in plots
Costep = 300; % Maximum number of iterations
nplots = 50; % Number of snapshots (plots) to take
plot_step = Costep/nplots; % Number of time steps between plots
dgap = 0.001; % Newton-Raphson interval
%--------------------------------------------------------------------------
%Preallocation and initial conditions
gap1_w =.5*ones(1,Nx); gap2_w =.5*ones(1,Nx);
u1_w = ones(1,Nx); u2_w = ones(1,Nx);
p1_w = zeros(1,Nx); gap1_wb = gap1_w;
%Setting up auxillary variables
for i=1:Nx
im(i)=i-1;
end
im(1)=1;
%--------------------------------------------------------------------------
for nt=1:Nt+1
t=(nt-1)*dt;
gap1_w(1) = .5*(1+.5*sin(pi*t)); % LHS BC for gap1_w
gap2_w(1) = .5*(1-.5*sin(pi*t)); % LHS BC for gap2_w
% From the mass conservation H1.u1 +H2.u2 =1 --------------------------
beta = (2+ cos(pi*t));
delta = (2+ cos(pi*t));
B = (gap1_w(1)*beta)+(gap2_w(1)*delta);
epsilon = 1/B;
%----------------------------------------------------------------------
u1_w(1) = epsilon*(2+ cos(pi*t)); % LHS BC for u1_w
u2_w(1) = epsilon*(2+ cos(pi*t)); % LHS BC for u2_w
p1_w(1) = epsilon*pi*(sin(pi*t)); % LHS BC for p1_w
%----------------------------------------------------------------------
% Save the values from the previous time-step as they are needed for the time derivative finite difference-
gap1n=gap1_w;
gap2n=gap2_w;
u1n=u1_w;
u2n=u2_w;
p1n=p1_w;
% Newton - Raphson for gap1_w + gap2_w = 1-----------------------------
flag=0;
mmm=1;
m=1;
gap1_w = gap1_wb;
gap2_w = 1-gap1_w;
while (flag == 0)
gap2_w = 1-gap1_w;
% Compute new gap1_w, gap2_w, u1_w, u2_w, p1_w using first order backward time finite difference scheme.
for i=1:Nx
% For gap 1----------------------------------------------------
gap1_w(i) = gap1n(i) - r*(u1_w(i).*(gap1_w(i)-gap1_w(im(i)))+ ...
gap1_w(i).*(u1_w(i)-u1_w(im(i))));
u1_w(i)=u1n(i)-r*(u1_w(i).*(u1_w(i)-u1_w(im(i)))+ (p1_w(i) - p1_w(im(i))));
% For gap 2----------------------------------------------------
gap2_w(i) = gap2n(i) - r*(u2_w(i).*(gap2_w(i)-gap2_w(im(i)))+ ...
gap2_w(i).*(u2_w(i)-u2_w(im(i))));
u2_w(i)=u2n(i)-r*(u2_w(i).*(u2_w(i)-u2_w(im(i))+(p1_w(i) - p1_w(im(i)))));
% Above, we have used p1 instead of p2 as they are equal.
end
for i=1:Nx
st(m,i) = 1-(gap1_w(i)+gap2_w(i));
end
m = m + 1;
if (m == 2)
gap1_w = gap1_w + dgap;
gap2_w = 1- gap1_w;
st(m)
end
if (m == 3)
for i=1:Nx
gap1_w(i) = (gap1_w(i).*st(1,i)-(gap1_w(i)-dgap).*st(2,i))./(st(1,i)-st(2,i));
end
gap2_w = 1- gap1_w;
end
if (m == 4)
mmm = mmm + 1;
if (mmm == 2)
m = 1;
else
flag = 1;
end
end
end
% Periodically record gap1_w, gap2_w, u1_w, u2_w, p1_w for ploting.
gap1_wplot(:,iplot) = gap1_w(:); % Record gap1_w(i) for plotting
gap2_wplot(:,iplot) = gap2_w(:); % Record gap2_w(i) for plotting
total(:,iplot) = gap1_w(:) + gap2_w(:); % Total gap for plotting
mass(:,iplot) = gap1_w(:).*u1_w(:) + gap2_w(:).*u2_w(:); % Total gap for plotting
u1_wplot(:,iplot) = u1_w(:); % Record u1_w(i) for plotting
u2_wplot(:,iplot) = u2_w(:); % Record u2_w(i) for plotting
p1_wplot(:,iplot) = p1_w(:); % Record p1_w(i) for plotting
tplot(iplot) = nt*dt; % Record time for plots
% end
iplot = iplot+1;
end
The thing is, it should use gap1_w,gap2_w u1_w,u2_w, p1_w on the right hand side in these lines using CURRENT TIME STEP NOT PREVIOUS
gap1_w(i) = gap1n(i) - r*(u1_w(i).*(gap1_w(i)-gap1_w(im(i)))+ ...
gap1_w(i).*(u1_w(i)-u1_w(im(i))));
u1_w(i)=u1n(i)-r*(u1_w(i).*(u1_w(i)-u1_w(im(i)))+ (p1_w(i) - p1_w(im(i))));
% For gap 2----------------------------------------------------
gap2_w(i) = gap2n(i) - r*(u2_w(i).*(gap2_w(i)-gap2_w(im(i)))+ ...
gap2_w(i).*(u2_w(i)-u2_w(im(i))));
u2_w(i)=u2n(i)-r*(u2_w(i).*(u2_w(i)-u2_w(im(i))+(p1_w(i) - p1_w(im(i)))
However, I assume it uses these from the previous time step. For example, instead of the gap1_w from the CURRENT time step, it recognises gap1_w as gap1n since it did not complete the loop.
Any help is much appreciated!

Answers (0)

Community Treasure Hunt

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

Start Hunting!