Solving a PDE with two variables using cank nicolson method

7 views (last 30 days)
Hana Bachi
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

Answers (2)

Alan Stevens
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 Comments
Torsten
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.

Sign in to comment.


Hana Bachi
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
  4 Comments
Hana Bachi
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

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!