# Solving a PDE with two variables using cank nicolson method

7 views (last 30 days)
Hana Bachi on 6 Apr 2022
Commented: Hana Bachi on 8 Apr 2022
Hello,
I'm trying to write a code using Crank Nicolson method to solve PDE with two varaiables. the problem is attached in the pdf file, problem N°:ii. Runing this code I'm getting an error message says:
"Array indices must be positive integers or logical values.
Error in Cp52nd (line 16)
f(x,y,t)= (x(i)+1)/((2*(1+y(j))*((1+t(k))^2)));"
the code is:
clear variables; close all;
%% Time and Space Discretization
%1. Space discretization
Lx = 1.0; dx =1/20; M=Lx/dx+1; x=0:dx:Lx;
Ly = 1.0; dy =1/20; N=Ly/dy+1; y=0:dy:Ly;
% Time discretization
tf =0.2; dt =1/20; O=tf/dt+1; t=0:dt:tf;
dx=0.1;
dy=0.1;
%Constants
for i=0:N
for j=0:M
for K=0:1
f(x,y,t)= (x(i)+1)/((2*(1+y(j))*((1+t(k))^2)));
g(x,y,t)= (y(j)+1)/((2*(1+x(i))*((1+t(k))^2)));
end
end
end
r=dt/(2*dx^2); a=2*r*f(x,y,t)+2*r*g(x,y,t)+1;c=-r*f(x,y,t);
U=zeros(M,N,O);
%% Loop to build the penthadiagonal A Matrix from Crank Nicolson Method
A=zeros((M-2)*(N-2),(M-2)*(N-2));
for j=1:(M-2)*(N-2)
for i=1:(M-2)*(N-2)
if i==j
A(i,j)=a;
elseif i-1==j && rem(i,N-2)~=1
A(i,j)=c;
elseif i+1==j && rem(i,N-2)~=0
A(i,j)=c;
elseif j-i==M-2
A(i,j)=c;
elseif i-j==N-2
A(i,j)=c;
end
end
end
%% Initial and boundary conditions from the analytical solution
%Initial conditions
for i=1:M
for j=1:N
U(i,j,1)=exp((x(i)+1)*(y(j)+1));
end
end
% Boundary conditions x=0 and x=1 for all y and t
for k=1:O
for j=1:N
U(1,j,k)=exp((x(1)+1)*(y(j)+1)*(t(k)+1));
U(M,j,k)=exp((x(M)+1)*(y(j)+1)*(t(k)+1));
end
end
% Boundary conditions y=0 and y=1 for all x and t
for k=1:O
for i=1:M
U(i,1,k)=exp((x(i)+1)*(y(1)+1)*(t(k)+1));
U(i,N,k)=exp((x(i)+1)*(y(N)+1)*(t(k)+1));
end
end
%% Solving the system Au=d
d=zeros((M-2)*(N-2),1);
%Building vector d at level n
for k=1:O-1
p=1;
for j=2:N-1
for i=2:M-1
d(p)=r*f(x(i),y(j),t(k))*U(i-1,j,k)+(1-2*r*f(x(i),y(j),t(k))-2*r*g(x(i),y(j),t(k)))*U(i,j,k)+r*f(x(i),y(j),t(k))*U(i+1,j,k)+r*g(x(i),y(j),t(k))*U(i,j-1,k)+r*g(x(i),y(j),t(k))*U(i,j+1,k);
if i==2
d(p)=d(p)+r*U(i-1,j,k+1);
elseif i==M-1
d(p)=d(p)+r*f(x(i),y(j),t(k))*U(i+1,j,k+1);
end
if j==2
d(p)=d(p)+r*g(x(i),y(j),t(k))*U(i,j-1,k+1);
elseif j==N-1
d(p)=d(p)+r*g(x(i),y(j),t(k))*U(i,j+1,k+1);
end
p=p+1;
end
end
%Solving u at level n+1
[u] = lu_dcmp(A,d); %Call function LU
q=1;
for j=2:N-1
for i=2:M-1
U(i,j,k+1)=u(q,1);
q=q+1;
end
end
end
%% Analytical solution and Absolut Error
for i=1:N
for j=1:M
UR(i,j)=exp((1+x(i))*(y(i)+1)*(1+tf));
Error(i,j)=abs(UR(i,j)-U(i,j,O)); %Absolut error
end
end
%% Plots
subplot(1,3,1)
mesh(x,y,U(:,:,O))
title('Numerical solution at t=0.2')
xlabel('x'); ylabel('y'); zlabel('U(x,y,t)');
%zlim([1 3])
subplot(1,3,2)
mesh(x,y,UR)
title('Analytical solution at t=0.2')
xlabel('x'); ylabel('y'); zlabel('U(x,y,t)');
%zlim([1 3])
subplot(1,3,3)
mesh(x,y,Error)
title('Error at t=0.2')
xlabel('x'); ylabel('y'); zlabel('Absolute Error');
%%
function [X] = lu_dcmp(A,b)
N=size(A,1);
p=bandwidth(A,'lower');
q=bandwidth(A,'upper');
L=zeros(N,N);
U=zeros(N,N);
for k=1:N-1
for i=k+1:min(k+p,N)
A(i,k)=A(i,k)/A(k,k);
for j=k+1:min(k+q,N)
A(i,j)=A(i,j)-A(i,k)*A(k,j);
end
end
end
for i=1:N
for j=1:N
if i>j
L(i,j)=A(i,j); %Creating the Lower Matrix
else
U(i,j)=A(i,j); %Creating the Upper Matrix
end
end
end
for i=1:N
sum=0;
for j=max(1,i-p):i-1
F=L(i,j)*b(j);
sum=sum+F;
end
b(i)=b(i)-sum;
end
X=zeros(N,1);
for i=1:N
sum=0;
i=N+1-i;
for j=i+1:min(i+q,N)
F=U(i,j)*X(j);
sum=sum+F;
end
X(i)=(b(i)-sum)/U(i,i);
end
end

Alan Stevens on 6 Apr 2022
Edited: Alan Stevens on 6 Apr 2022
Matlab indices start at 1, so the loops i, j and K in the following
for i=0:N
for j=0:M
for K=0:1
f(x,y,t)= (x(i)+1)/((2*(1+y(j))*((1+t(k))^2)));
g(x,y,t)= (y(j)+1)/((2*(1+x(i))*((1+t(k))^2)));
end
end
end
are invalid when i, j K equal 0.
Also, you probably need
f(i,j,K) = ...
etc.
And be consistent with the case: do you want K or k?
##### 2 CommentsShowHide 1 older comment
Torsten on 6 Apr 2022
Shift all your arrays one position to the right to start at i,j,k = 1.
Then your boundary condition is in position 1 instead of 0 and your end point in position n+1 instead of n - it doesn't matter.
The error message is the consequence that your loops start aat 0 instead of 1.

Hana Bachi on 7 Apr 2022
I modified the arrays and relpaced f(x,y,t) and g(x,y,t) by their expresions. I'm getting a plot but the analytical solution looks totally different from the numerical solution. where is the mistake with this code:
the plot is in below
the enw code:
%%Numerical Methods CP5
%Problem 2
clear variables; close all;
%% Time and Space Discretization
%1. Space discretization
Lx = 1.0; dx =1/20; M=Lx/dx+1; x=0:dx:Lx;
Ly = 1.0; dy =1/20; N=Ly/dy+1; y=0:dy:Ly;
% Time discretization
tf =0.2; dt =1/20; O=tf/dt+1; t=0:dt:tf;
dx=0.1;
dy=0.1;
%Constants
for i=1:N
for j=1:M
for k=1:1+1
f(i,j,k)= (x(i)+1)/((2*(1+y(j))*((1+t(k))^2)));
g(i,j,k)= (y(j)+1)/((2*(1+x(i))*((1+t(k))^2)));
end
end
end
r=dt/(2*dx^2); a=2*r*(x(i)+1)/((2*(1+y(j))*((1+t(k))^2)))+2*r*(y(j)+1)/((2*(1+x(i))*((1+t(k))^2)))+1;c=-r*(x(i)+1)/((2*(1+y(j))*((1+t(k))^2)));
U=zeros(M,N,O);
%% Loop to build the penthadiagonal A Matrix from Crank Nicolson Method
A=zeros((M-2)*(N-2),(M-2)*(N-2));
for j=1:(M-2)*(N-2)
for i=1:(M-2)*(N-2)
if i==j
A(i,j)=a;
elseif i-1==j && rem(i,N-2)~=1
A(i,j)=c;
elseif i+1==j && rem(i,N-2)~=0
A(i,j)=c;
elseif j-i==M-2
A(i,j)=c;
elseif i-j==N-2
A(i,j)=c;
end
end
end
%% Initial and boundary conditions from the analytical solution
%Initial conditions
for i=1:M
for j=1:N
U(i,j,1)=exp((x(i)+1)*(y(j)+1));
end
end
% Boundary conditions x=0 and x=1 for all y and t
for k=1:O
for j=1:N
U(1,j,k)=exp((x(1)+1)*(y(j)+1)*(t(k)+1));
U(M,j,k)=exp((x(M)+1)*(y(j)+1)*(t(k)+1));
end
end
% Boundary conditions y=0 and y=1 for all x and t
for k=1:O
for i=1:M
U(i,1,k)=exp((x(i)+1)*(y(1)+1)*(t(k)+1));
U(i,N,k)=exp((x(i)+1)*(y(N)+1)*(t(k)+1));
end
end
%% Solving the system Au=d
d=zeros((M-2)*(N-2),1);
%Building vector d at level n
for k=1:O-1
p=1;
for j=2:N-1
for i=2:M-1
d(p)=r*(x(i)+1)/((2*(1+y(j))*((1+t(k))^2)))*U(i-1,j,k)+(1-2*r*(x(i)+1)/((2*(1+y(j))*((1+t(k))^2))))-2*r*(y(j)+1)/((2*(1+x(i))*((1+t(k))^2)))*U(i,j,k)+r*(x(i)+1)/((2*(1+y(j))*((1+t(k))^2)))*U(i+1,j,k)+r*(y(j)+1)/((2*(1+x(i))*((1+t(k))^2)))+r*(y(j)+1)/((2*(1+x(i))*((1+t(k))^2)))*U(i,j+1,k);
if i==2
d(p)=d(p)+r*U(i-1,j,k+1);
elseif i==M-1
d(p)=d(p)+r*(x(i)+1)/((2*(1+y(j))*((1+t(k))^2)))*U(i+1,j,k+1);
end
if j==2
d(p)=d(p)+r*(y(j)+1)/((2*(1+x(i))*((1+t(k))^2)))*U(i,j-1,k+1);
elseif j==N-1
d(p)=d(p)+r*(y(j)+1)/((2*(1+x(i))*((1+t(k))^2)))*U(i,j+1,k+1);
end
p=p+1;
end
end
%Solving u at level n+1
[u] = lu_dcmp(A,d); %Call function LU
q=1;
for j=2:N-1
for i=2:M-1
U(i,j,k+1)=u(q,1);
q=q+1;
end
end
end
%% Analytical solution and Absolut Error
for i=1:N
for j=1:M
UR(i,j)=exp((1+x(i))*(y(i)+1)*(1+tf));
Error(i,j)=abs(UR(i,j)-U(i,j,O)); %Absolut error
end
end
%% Plots
subplot(1,3,1)
mesh(x,y,U(:,:,O))
title('Numerical solution at t=0.2')
xlabel('x'); ylabel('y'); zlabel('U(x,y,t)');
%zlim([1 3])
subplot(1,3,2)
mesh(x,y,UR)
title('Analytical solution at t=0.2')
xlabel('x'); ylabel('y'); zlabel('U(x,y,t)');
%zlim([1 3])
subplot(1,3,3)
mesh(x,y,Error)
title('Error at t=0.2')
xlabel('x'); ylabel('y'); zlabel('Absolute Error');
%%
function [X] = lu_dcmp(A,b)
N=size(A,1);
p=bandwidth(A,'lower');
q=bandwidth(A,'upper');
L=zeros(N,N);
U=zeros(N,N);
for k=1:N-1
for i=k+1:min(k+p,N)
A(i,k)=A(i,k)/A(k,k);
for j=k+1:min(k+q,N)
A(i,j)=A(i,j)-A(i,k)*A(k,j);
end
end
end
for i=1:N
for j=1:N
if i>j
L(i,j)=A(i,j); %Creating the Lower Matrix
else
U(i,j)=A(i,j); %Creating the Upper Matrix
end
end
end
for i=1:N
sum=0;
for j=max(1,i-p):i-1
F=L(i,j)*b(j);
sum=sum+F;
end
b(i)=b(i)-sum;
end
X=zeros(N,1);
for i=1:N
sum=0;
i=N+1-i;
for j=i+1:min(i+q,N)
F=U(i,j)*X(j);
sum=sum+F;
end
X(i)=(b(i)-sum)/U(i,i);
end
end
Hana Bachi on 8 Apr 2022
Yes, I used crank nicolson descritization but I'm still getting different results. I'm not sure if the mistake is with the BC, IC or d(p) itself. Here is my last results