Solving a PDE with two variables using cank nicolson method

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

Answer 1

[Alan Stevens]

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?

[Hana Bachi]

1- how can I define i, j and k at 0 because it's boundary condition.

2- After changing f(x,y,t) and g(x,y,t) to f(i,j,K) and g(i,j,K), I got this errro

Array indices must be positive integers or logical values.

Error in Cp52nd (line 16)

f(i,j,k)= (x(i)+1)/((2*(1+y(j))*((1+t(k))^2)));

[Torsten]

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.

Answer 2

[Hana Bachi]

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

[Torsten]

Debugging - the hardest time ...

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)));

This won't work since it must be done within a loop over i and j - in the loop where you define the matrix elements.

I suggest you write out the equations for Crank-Nicolson first on paper.

Then you can see how you have to define the matrices and vectors of the right-hand side.

Start with

(u_ij^(k+1) - u_ij^(k)) / dt = 1/2 * (f(x_i,y_j,t_k,u_ij^(k)) + f(x_i,y_j,t_(k+1),u_ij^(k+1))

where

f(x_i,y_j,t_k,u_ij) = (x_i+1)/(2*(1+y_j)*(1+t_k)^2 ) * (u_(i-1),j - 2*u_ij + u_(i+1),j) / dx^2 +
                      (1+y_j)/(2*(1+x_i)*(1+t_k)^2 ) * (u_i,(j-1) - 2*u_ij + u_i,(j+1)) / dy^2     

[Hana Bachi]

I don't understand how you write f(x,y,t,u) while U is the one which is on the function of f,g,x,y,t ( U(x,y,t,f,g))

[Torsten]

I didn't use variable names from your code.

I just wrote down which equations you have to implement for Crank-Nicolson.

In inner grid points, these are

(u_ij^(k+1) - u_ij^(k)) / dt = 1/2 * (f(x_i,y_j,t_k,u_ij^(k)) + f(x_i,y_j,t_(k+1),u_ij^(k+1))

with

f(x_i,y_j,t_k,u_ij) = (x_i+1)/(2*(1+y_j)*(1+t_k)^2 ) * (u_(i-1),j - 2*u_ij + u_(i+1),j) / dx^2 +
                      (1+y_j)/(2*(1+x_i)*(1+t_k)^2 ) * (u_i,(j-1) - 2*u_ij + u_i,(j+1)) / dy^2
                      

In boundary points, these are (e.g. for x=0):

 0.5*(u_1j^(k) + u_1j^(k+1)) = 0.5*(exp((1+y_j)*(1+t_k))+exp((1+y_j)*(1+t_(k+1))))

[Hana Bachi]

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