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

6 views (last 30 days)
Meva on 21 Feb 2016
Edited: Meva on 21 Feb 2016
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!