Solving a PDE with two variables using cank nicolson method

2 visualizaciones (últimos 30 días)
Hana Bachi
Hana Bachi el 6 de Abr. de 2022
Respondida: Kuldeep Malik el 28 de Jul. de 2023
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

Respuestas (3)

Alan Stevens
Alan Stevens el 6 de Abr. de 2022
Editada: Alan Stevens el 6 de Abr. de 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 comentarios
Hana Bachi
Hana Bachi el 6 de Abr. de 2022
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
Torsten el 6 de Abr. de 2022
Editada: Torsten el 6 de Abr. de 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.

Iniciar sesión para comentar.


Hana Bachi
Hana Bachi el 7 de Abr. de 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 comentarios
Torsten
Torsten el 7 de Abr. de 2022
Editada: Torsten el 7 de Abr. de 2022
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
Hana Bachi el 8 de Abr. de 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

Iniciar sesión para comentar.


Kuldeep Malik
Kuldeep Malik el 28 de Jul. de 2023
I am stuck with the code of solving pde using Crank nicolson method although code is running well but numerical approach is way difference than exact solution
is something wrong with code?/
Help Please
%Numerical solution except Initial and Boundary Values
clc, clear, close
tic
% Parameters
% lambda
dx=0.1;% discretization size for x axis
dt=0.01;
a=0.5;
L=1;
La=(dx/dt)^a;
nx=11; % or nx=[(10-0)/dx]+1 with L=1
n=3;
for nt=2:n
% to compute 10 time steps k=[1,11]
A=zeros(9);
%Define Initial and Boundary conditions
for j=1:nt
u(1,j)=-((mlf(a,1,((j-1).*dt).^a,7)-mlf(a,1,-1.*(((j-1).*dt).^a),7))./2);
u(nx,j)=((mlf(a,1,((nx-1).*dx).^a,7)-mlf(a,1,-1.*(((nx-1).*dx).^a),7))./2).*((mlf(a,1,((j-1).*dt).^a,7)+mlf(a,1,-1.*((j-1).*dt).^a,7))./2)-((mlf(a,1,((nx-1).*dx).^a,7)+mlf(a,1,-1.*(((nx-1).*dx).^a),7))./2).*((mlf(a,1,((j-1).*dt).^a,7)-mlf(a,1,-1.*((j-1).*dt).^a,7))./2);
end
for i=2:nx-1
u(i,1)=(mlf(a,1,((i-1).*dx).^a,7)-mlf(a,1,-1.*(((i-1).*dx).^a),7))./2;
end
%Defining omega values
for ii=1:nx-1
wa(ii)=(((-1).^ii).*gamma(a+1))./(gamma(ii+1).*gamma(a-ii+1));
end
wa;
%Defining A matrix
A=zeros(nx-2);
for m=1:nx-2
for mm=1:nx-2
if(m==mm)
A(m,mm)=La*wa(1)+1;
elseif(m==mm+1)
A(m,mm)=wa(2);
elseif(m==mm+2)
A(m,mm)=wa(3);
elseif(m==mm+3)
A(m,mm)=wa(4);
elseif(m==mm+4)
A(m,mm)=wa(5);
elseif(m==mm+5)
A(m,mm)=wa(6);
elseif(m==mm+6)
A(m,mm)=wa(7);
elseif(m==mm+7)
A(m,mm)=wa(8);
elseif(m==mm+8)
A(m,mm)=wa(9);
else
A(m,mm)=0;
end
end
end
A;
%Defining B matrix
summ=zeros(9,1);
for pp=2:nx-1
sum=0;
for kk=1:nt
sum=sum+wa(kk+1)*u(pp,nt+1-kk);
end
summ(pp,1)=sum;
end
summ;
B=zeros(nx-2,1);
for jj=1:nx-2
B(jj,1)=-La*summ(jj,1)-wa(jj+1)*u(1,nt);
end
B;
U(nt,:)=(A\B)';
end
Num=U
%Exact solution
for jjj=1:nt
for iii=1:nx
% u(1,j)=-(mlf(a,1,((j-1).*dt).^a,a)-mlf(a,1,-1.*((j-1).*dt).^a,a))./2;
Ex(iii,jjj)=((mlf(a,1,((iii-1).*dx).^a,7)-mlf(a,1,-1.*(((iii-1).*dx).^a),7))./2).*((mlf(a,1,((jjj-1).*dt).^a,7)+mlf(a,1,-1.*((jjj-1).*dt).^a,7))./2)-((mlf(a,1,((iii-1).*dx).^a,7)+mlf(a,1,-1.*(((iii-1).*dx).^a),7))./2).*((mlf(a,1,((jjj-1).*dt).^a,7)-mlf(a,1,-1.*((jjj-1).*dt).^a,7))./2);
end
end
Exact=Ex'
All the related t terms are there

Productos

Community Treasure Hunt

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

Start Hunting!

Translated by