Solving a PDE with two variables using cank nicolson method
3 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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
0 commentaires
Réponses (3)
Alan Stevens
le 6 Avr 2022
Modifié(e) : Alan Stevens
le 6 Avr 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 commentaires
Torsten
le 6 Avr 2022
Modifié(e) : Torsten
le 6 Avr 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
le 7 Avr 2022
4 commentaires
Torsten
le 7 Avr 2022
Modifié(e) : Torsten
le 7 Avr 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))))
Kuldeep Malik
le 28 Juil 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'
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1444282/image.jpeg)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1444287/image.jpeg)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1444292/image.jpeg)
All the related t terms are there
0 commentaires
Voir également
Produits
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!