THANK YOU!!!! For those wondering about dx ~= dy line 25 == -2*( (dx/dy)^2 +1) line 28&29== (dx/dy)^2 I believe this is all that is needed. I could be wrong. I'm a rookie.
Sir, this video is really helpful. I can understand everything. Sir, Please make videos for parabolic and hyperbolic equations also. And I have one doubt. Why we needed transpose the matrix while plotting??
great explanation God bless you....i have been looking for such video.if you are volenteer can i get your program and tray it with other boundry condition
what should one do if there are finite boundary values instead of x or 1-y? Does that mean I can not create a grid with size greater that the number of given boundary values? Thank you for the code. It was really amazing
I realized that the equation is not elliptic but parabolic with 2D (spatial) and time dependent. The paper which I am implementing says that they found the initial values using a 5-degree polynomial (for which they have shown the values of initial condition (and BC) in a table for 7x7 grid but in results they have used a grid size of 1000x1000 and the polynomial equation does not satisfy the initial values (or BC) in their table when I plug in the x and y in f(x,y) (which is a given 5 degree polynomial)
Hi... i just wanted to ask how would the coding change if i was using fourth order differential equation to solve for a plate with variable support conditions (boundary)
I totally appreciate your work and I'm trying to do exactly the same code even with the same symbols but it doesn't work and it says" Warning: Matrix is singular to working precision" Do you have an idea what's wrong. Thanks : )
+Hassan Saleem This error occurs when there is a row in the M matrix that has all zeros. If you make your grid very coarse, like 4x4, and print the M matrix you will likely find a row of all zeros. Then you can figure out which unknown this row corresponds to and track down the error. Most likely one of the boundary conditions is wrong.
@@AhmedAbdelazim-- Nx=50; Ny=50; Lx=5; Ly=5; x=linspace(0,Lx,Nx); y=linspace(0,Ly,Ny); dx=x(2)-x(1); dy=y(2)-y(1); N=Nx*Ny; M=zeros(N,N); B=zeros(N,1); %phi_vec=inv(M); for i=2:Nx-1 for j=2:Ny-1 n=i+(j-1)*Nx; M(n,n)= -4; M(n,n-1)= 1; M(n,n+1)= 1; M(n,n-Nx)=1; M(n,n+Nx)=1; B(n,1) =0; end end % BC: i=1; for j=1:Ny; n=i+(j-1)*Nx;%nth linha para essa ij M(n,n)=1;%diagonal principal B(n,1)=y(j);%CI am y_j end %lado direito phi=1-y i=Nx; for j=1:Ny n=i+(j-1)*Nx;%nth linha para essa ij M(n,n)=1;%diagonal principal B(n,1)=1-y(j);%CI am y_j end %lado de baixo phi=x j=1; for i=1:Nx n=i+(j-1)*Nx;%nth linha para essa ij M(n,n)=1;%diagonal principal B(n,1)=x(i);%CI am x_i end %lado de cima phi=1-x j=Ny; for i=1:Nx n=i+(j-1)*Nx;%nth linha para essa ij M(n,n)=1;%diagonal principal B(n,1)=1-x(i);%CI am x_i end %phi_vec=M\B; disp(M); phi_vec=inv(M)*B; for i=1:Nx for j=1:Ny n=i+(j-1)*Nx; phi(i,j)=phi_vec(n); end end surf(x,y,phi') xlabel('x') ylabel('y')
i wrote the same at Matlab2020 but this error becomes: "Error using / Matrix dimensions must agree. Error in EllipticEq (line 62) phi_vec = M/B;" can you help?
It is the conversion of 2D vector into 1D vector. Like this, A = 11 12 22 23 dimension of (A) = 2,2 i.e i = 1,2 and j = 1,2 to get A(2,2) = 23, when it is written in a 1D vector as -- vec(A) = 11, 12, 22, 23 n = 2 + (2-1)*2 = 2 + 1*2 = 4th element in vec(A) = 23. Hope it helps someone. :)
i write this code and my graphic inst like your what i am doing wrong. %Elliptic PDE-FiniteDifference 2D clear clc %Imputs Nx=50; Ny=50; Lx=1; Ly=1; %grid x=linspace(0,Lx,Nx); y=linspace(0,Ly,Ny); dx=x(2)-x(1);%mesh size dy=y(2)-y(1); %initialize matrices N=Nx*Ny;%numeros de variaveis M=zeros(N,N);%N linhas, N colunas B=zeros(N,1);%N linhas, 1 coluna %Pontos interiores for i=2:Nx-1% loop para direçao x , nao levando em consideraçao o primeiro for j=2:Ny-1% loop para direçao y , nao levando em consideraçao o primeiro n=i+(j-1)*Nx; %convert ij grid point to the nth M(n,n)=-4;%diagonal principal M(n,n-1)=1; M(n,n+1)=1; M(n,n-Nx)=1; M(n,n+Nx)=1; B(n,1)=0; end end %Condiçoes de contorno %lado esquerdo phi=y i=1; for j=1:Ny; n=i+(j-1)*Nx;%nth linha para essa ij M(n,n)=1;%diagonal principal B(n,1)=y(j);%CI am y_j end %lado direito phi=1-y i=Nx; for j=1:Ny n=i+(j-1)*Nx;%nth linha para essa ij M(n,n)=1;%diagonal principal B(n,1)=1-y(j);%CI am y_j end %lado de baixo phi=x j=1; for i=1:Nx n=i+(j-1)*Nx;%nth linha para essa ij M(n,n)=1;%diagonal principal B(n,1)=x(i);%CI am x_i end %lado de cima phi=1-x j=Ny; for i=1:Nx n=i+(j-1)*Nx;%nth linha para essa ij M(n,n)=1;%diagonal principal B(n,1)=1-x(i);%CI am x_i end %Resolvendo phi_vec=M\B; %Convertendo vetor para uma resposta 2d for i=1:Nx for j=1:Ny n=1+(j-1)*Nx;%nth linha para essa ij phi(i,j)=phi_vec(n); end end %plotando resultado surf(x,y,phi')%' transposta xlabel('x') ylabel('y') set(gca,'Fontsize',16) %computando velocidade do campo potencial %u=-dphi/dx% componente x da velocidade %v=-dphi/dy% componente y da velocidade u=zeros(Nx,Ny); v=zeros(Nx,Ny); for i=2:Nx-1% loop para os pontos interiores for j=2:Ny-1 u(i,j)=-(phi(i+1,j)-phi(i-1,j))/(2*dx); v(i,j)=-(phi(i,j+1)-phi(i,j-1))/(2*dy); end end figure(2);clf(2) quiver(x,y,u',v') xlabel('x') ylabel('y') title('Campo de velocidade')
fixed it: Nx=50; Ny=50; Lx=5; Ly=5; x=linspace(0,Lx,Nx); y=linspace(0,Ly,Ny); dx=x(2)-x(1); dy=y(2)-y(1); N=Nx*Ny; M=zeros(N,N); B=zeros(N,1); %phi_vec=inv(M); for i=2:Nx-1 for j=2:Ny-1 n=i+(j-1)*Nx; M(n,n)= -4; M(n,n-1)= 1; M(n,n+1)= 1; M(n,n-Nx)=1; M(n,n+Nx)=1; B(n,1) =0; end end % BC: i=1; for j=1:Ny; n=i+(j-1)*Nx;%nth linha para essa ij M(n,n)=1;%diagonal principal B(n,1)=y(j);%CI am y_j end %lado direito phi=1-y i=Nx; for j=1:Ny n=i+(j-1)*Nx;%nth linha para essa ij M(n,n)=1;%diagonal principal B(n,1)=1-y(j);%CI am y_j end %lado de baixo phi=x j=1; for i=1:Nx n=i+(j-1)*Nx;%nth linha para essa ij M(n,n)=1;%diagonal principal B(n,1)=x(i);%CI am x_i end %lado de cima phi=1-x j=Ny; for i=1:Nx n=i+(j-1)*Nx;%nth linha para essa ij M(n,n)=1;%diagonal principal B(n,1)=1-x(i);%CI am x_i end %phi_vec=M\B; disp(M); phi_vec=inv(M)*B; for i=1:Nx for j=1:Ny n=i+(j-1)*Nx; phi(i,j)=phi_vec(n); end end surf(x,y,phi') xlabel('x') ylabel('y')
I have been looking for this kind of explanation everywhere. I finally found it here. Thank you so much.
This is the most concise explanation of the discretization and the solution of a PDE on youtube, thank you
Wow Great video
It was very helpful for my project sir. Thanks a lot.
THANK YOU!!!!
For those wondering about dx ~= dy
line 25 == -2*( (dx/dy)^2 +1)
line 28&29== (dx/dy)^2
I believe this is all that is needed. I could be wrong. I'm a rookie.
Thank you sir, your explanation regarding the Matlab code cleared up several issues I had
nice videos. this is such a great help that i can understand matlab PDE coding.
I would like to thank you so so much ! you don't imagine how hard you saved my life
I used it in my Heat transfer course project :) I faced lots of difficulties, but thankfully I could solve them
@@ath216 can you helpe me
@@lobnalobna7489 what kind of help do you need, I almost forgot everything about this topic, But I can try..
و انا ايضا احب الرياضيات :)
Sir You are really great
Sir, this video is really helpful. I can understand everything. Sir, Please make videos for parabolic and hyperbolic equations also.
And I have one doubt. Why we needed transpose the matrix while plotting??
how if the boundary condition for top is dphi/dy=0? what should be the coding for this boundary condition?
I thoroughly enjoyed listening to and learning from this comprehensive lesson, keep up the great work are there more of this?
great explanation God bless you....i have been looking for such video.if you are volenteer can i get your program and tray it with other boundry condition
what should one do if there are finite boundary values instead of x or 1-y? Does that mean I can not create a grid with size greater that the number of given boundary values?
Thank you for the code. It was really amazing
If you have a finite value, something like phi(x,y=0)=10. Then, just replace phi=-x with phi=10. You can still have a fine grid.
I realized that the equation is not elliptic but parabolic with 2D (spatial) and time dependent. The paper which I am implementing says that they found the initial values using a 5-degree polynomial (for which they have shown the values of initial condition (and BC) in a table for 7x7 grid but in results they have used a grid size of 1000x1000 and the polynomial equation does not satisfy the initial values (or BC) in their table when I plug in the x and y in f(x,y) (which is a given 5 degree polynomial)
Sir, could you please do the same for parabolic PDEs anytime soon?
Sir, please solve the same using finite element method.
Thank you
is this a 2D elliptic equation?
any one who can help me to how to formulate the bounder condition if the flow region is L shaped rather than square
Hi... i just wanted to ask how would the coding change if i was using fourth order differential equation to solve for a plate with variable support conditions (boundary)
if you got any idea or help from anyone to solve four oder ode than send me
the reason why you got a wrong solution initially is because you calculated phi=M/B but should have been phi=B/M
I totally appreciate your work and I'm trying to do exactly the same code even with the same symbols but it doesn't work and it says" Warning: Matrix is singular to working precision" Do you have an idea what's wrong. Thanks : )
+Hassan Saleem This error occurs when there is a row in the M matrix that has all zeros. If you make your grid very coarse, like 4x4, and print the M matrix you will likely find a row of all zeros. Then you can figure out which unknown this row corresponds to and track down the error. Most likely one of the boundary conditions is wrong.
+Esraa Mohamed i a problem like you
i figured out what was wrong with my code and the code in the vedio is totally right...
@@esraamohamed-zy6cn can u please tell me what was your mistake ?? as the same is happenning with me now
@@AhmedAbdelazim-- Nx=50;
Ny=50;
Lx=5;
Ly=5;
x=linspace(0,Lx,Nx);
y=linspace(0,Ly,Ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
N=Nx*Ny;
M=zeros(N,N);
B=zeros(N,1);
%phi_vec=inv(M);
for i=2:Nx-1
for j=2:Ny-1
n=i+(j-1)*Nx;
M(n,n)= -4;
M(n,n-1)= 1;
M(n,n+1)= 1;
M(n,n-Nx)=1;
M(n,n+Nx)=1;
B(n,1) =0;
end
end
% BC:
i=1;
for j=1:Ny;
n=i+(j-1)*Nx;%nth linha para essa ij
M(n,n)=1;%diagonal principal
B(n,1)=y(j);%CI am y_j
end
%lado direito phi=1-y
i=Nx;
for j=1:Ny
n=i+(j-1)*Nx;%nth linha para essa ij
M(n,n)=1;%diagonal principal
B(n,1)=1-y(j);%CI am y_j
end
%lado de baixo phi=x
j=1;
for i=1:Nx
n=i+(j-1)*Nx;%nth linha para essa ij
M(n,n)=1;%diagonal principal
B(n,1)=x(i);%CI am x_i
end
%lado de cima phi=1-x
j=Ny;
for i=1:Nx
n=i+(j-1)*Nx;%nth linha para essa ij
M(n,n)=1;%diagonal principal
B(n,1)=1-x(i);%CI am x_i
end
%phi_vec=M\B;
disp(M);
phi_vec=inv(M)*B;
for i=1:Nx
for j=1:Ny
n=i+(j-1)*Nx;
phi(i,j)=phi_vec(n);
end
end
surf(x,y,phi')
xlabel('x')
ylabel('y')
i wrote the same at Matlab2020 but this error becomes: "Error using /
Matrix dimensions must agree.
Error in EllipticEq (line 62)
phi_vec = M/B;" can you help?
You need a back slash not a forward slash.
@@prof.markowkes-montanastat3851 i fixed it thank you sir
How can I define the boundary in my case 1000, 500, 50, and 0 ?. how can I apply it in the code ? some one please help me :)
please explain this a little n=i+( j-1) * nx
I don't understand only this, further more these are very very helpful videos....Thanks a lot
It is the conversion of 2D vector into 1D vector. Like this,
A = 11 12
22 23
dimension of (A) = 2,2 i.e i = 1,2 and j = 1,2
to get A(2,2) = 23, when it is written in a 1D vector as -- vec(A) = 11, 12, 22, 23
n = 2 + (2-1)*2 = 2 + 1*2 = 4th element in vec(A) = 23.
Hope it helps someone. :)
what will be formula for non uniform grid?
@@fadoobaba I believe the same. I modified the code and left that the same and I got almost exact solutions.
please send the link of pdf file of matlab code
i write this code and run i got singularity
error
i write this code and my graphic inst like your what i am doing wrong.
%Elliptic PDE-FiniteDifference 2D
clear
clc
%Imputs
Nx=50;
Ny=50;
Lx=1;
Ly=1;
%grid
x=linspace(0,Lx,Nx);
y=linspace(0,Ly,Ny);
dx=x(2)-x(1);%mesh size
dy=y(2)-y(1);
%initialize matrices
N=Nx*Ny;%numeros de variaveis
M=zeros(N,N);%N linhas, N colunas
B=zeros(N,1);%N linhas, 1 coluna
%Pontos interiores
for i=2:Nx-1% loop para direçao x , nao levando em consideraçao o primeiro
for j=2:Ny-1% loop para direçao y , nao levando em consideraçao o primeiro
n=i+(j-1)*Nx; %convert ij grid point to the nth
M(n,n)=-4;%diagonal principal
M(n,n-1)=1;
M(n,n+1)=1;
M(n,n-Nx)=1;
M(n,n+Nx)=1;
B(n,1)=0;
end
end
%Condiçoes de contorno
%lado esquerdo phi=y
i=1;
for j=1:Ny;
n=i+(j-1)*Nx;%nth linha para essa ij
M(n,n)=1;%diagonal principal
B(n,1)=y(j);%CI am y_j
end
%lado direito phi=1-y
i=Nx;
for j=1:Ny
n=i+(j-1)*Nx;%nth linha para essa ij
M(n,n)=1;%diagonal principal
B(n,1)=1-y(j);%CI am y_j
end
%lado de baixo phi=x
j=1;
for i=1:Nx
n=i+(j-1)*Nx;%nth linha para essa ij
M(n,n)=1;%diagonal principal
B(n,1)=x(i);%CI am x_i
end
%lado de cima phi=1-x
j=Ny;
for i=1:Nx
n=i+(j-1)*Nx;%nth linha para essa ij
M(n,n)=1;%diagonal principal
B(n,1)=1-x(i);%CI am x_i
end
%Resolvendo
phi_vec=M\B;
%Convertendo vetor para uma resposta 2d
for i=1:Nx
for j=1:Ny
n=1+(j-1)*Nx;%nth linha para essa ij
phi(i,j)=phi_vec(n);
end
end
%plotando resultado
surf(x,y,phi')%' transposta
xlabel('x')
ylabel('y')
set(gca,'Fontsize',16)
%computando velocidade do campo potencial
%u=-dphi/dx% componente x da velocidade
%v=-dphi/dy% componente y da velocidade
u=zeros(Nx,Ny);
v=zeros(Nx,Ny);
for i=2:Nx-1% loop para os pontos interiores
for j=2:Ny-1
u(i,j)=-(phi(i+1,j)-phi(i-1,j))/(2*dx);
v(i,j)=-(phi(i,j+1)-phi(i,j-1))/(2*dy);
end
end
figure(2);clf(2)
quiver(x,y,u',v')
xlabel('x')
ylabel('y')
title('Campo de velocidade')
fixed it:
Nx=50;
Ny=50;
Lx=5;
Ly=5;
x=linspace(0,Lx,Nx);
y=linspace(0,Ly,Ny);
dx=x(2)-x(1);
dy=y(2)-y(1);
N=Nx*Ny;
M=zeros(N,N);
B=zeros(N,1);
%phi_vec=inv(M);
for i=2:Nx-1
for j=2:Ny-1
n=i+(j-1)*Nx;
M(n,n)= -4;
M(n,n-1)= 1;
M(n,n+1)= 1;
M(n,n-Nx)=1;
M(n,n+Nx)=1;
B(n,1) =0;
end
end
% BC:
i=1;
for j=1:Ny;
n=i+(j-1)*Nx;%nth linha para essa ij
M(n,n)=1;%diagonal principal
B(n,1)=y(j);%CI am y_j
end
%lado direito phi=1-y
i=Nx;
for j=1:Ny
n=i+(j-1)*Nx;%nth linha para essa ij
M(n,n)=1;%diagonal principal
B(n,1)=1-y(j);%CI am y_j
end
%lado de baixo phi=x
j=1;
for i=1:Nx
n=i+(j-1)*Nx;%nth linha para essa ij
M(n,n)=1;%diagonal principal
B(n,1)=x(i);%CI am x_i
end
%lado de cima phi=1-x
j=Ny;
for i=1:Nx
n=i+(j-1)*Nx;%nth linha para essa ij
M(n,n)=1;%diagonal principal
B(n,1)=1-x(i);%CI am x_i
end
%phi_vec=M\B;
disp(M);
phi_vec=inv(M)*B;
for i=1:Nx
for j=1:Ny
n=i+(j-1)*Nx;
phi(i,j)=phi_vec(n);
end
end
surf(x,y,phi')
xlabel('x')
ylabel('y')