Elliptic PDE - FiniteDifference - Part 3 - MATLAB code

[@artcellCTRL]

This is the most concise explanation of the discretization and the solution of a PDE on youtube, thank you

[@w_wahlang1990]

I have been looking for this kind of explanation everywhere. I finally found it here. Thank you so much.

[@saidur_abir]

It was very helpful for my project sir. Thanks a lot.

[@ath216]

I would like to thank you so so much ! you don't imagine how hard you saved my life

[@ath216]

I used it in my Heat transfer course project :) I faced lots of difficulties, but thankfully I could solve them

[@lobnalobna7489]

@@ath216 can you helpe me

[@ath216]

@@lobnalobna7489 what kind of help do you need, I almost forgot everything about this topic, But I can try.. و انا ايضا احب الرياضيات :)

[@alexnguyennn]

nice videos. this is such a great help that i can understand matlab PDE coding.

[@nineteenmoscow3829]

Thank you sir, your explanation regarding the Matlab code cleared up several issues I had

[@walkerhendricks5754]

Wow Great video

[@ongyongshien5844]

how if the boundary condition for top is dphi/dy=0? what should be the coding for this boundary condition?

[@jasonbowens8369]

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.

[@rahulsinha7192]

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??

[@10010ist]

I thoroughly enjoyed listening to and learning from this comprehensive lesson, keep up the great work are there more of this?

[@yohanneskassa3786]

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

[@biplabhembram3042]

Sir You are really great

[@yohanneskassa3786]

any one who can help me to how to formulate the bounder condition if the flow region is L shaped rather than square

[@nihalrajendra6411]

the reason why you got a wrong solution initially is because you calculated phi=M/B but should have been phi=B/M

[@ehsanzahoor607]

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)

[@naseemullahkhan9008]

if you got any idea or help from anyone to solve four oder ode than send me

[@gohteikwei1289]

is this a 2D elliptic equation?

[@XclusiVvHD]

Sir, could you please do the same for parabolic PDEs anytime soon?

[@NithinGoona]

Sir, please solve the same using finite element method. Thank you

[@MuhammadArshad]

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

[@prof.markowkes-montanastat3851]

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.

[@MuhammadArshad]

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)

[@elifmath1229]

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?

[@prof.markowkes-montanastat3851]

You need a back slash not a forward slash.

[@elifmath1229]

@@prof.markowkes-montanastat3851 i fixed it thank you sir

[@esraamohamed-zy6cn]

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

[@prof.markowkes-montanastat3851]

+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.

[@oussamadraidi5275]

+Esraa Mohamed i a problem like you

[@esraamohamed-zy6cn]

i figured out what was wrong with my code and the code in the vedio is totally right...

[@AhmedAbdelazim--]

​@@esraamohamed-zy6cn can u please tell me what was your mistake ?? as the same is happenning with me now

[@muhammadusama6173]

@@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')

[@syaffelectricaldesigner6534]

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

[@AsadKhan-kg5yk]

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

[@MuhammadArshad]

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. :)

[@fadoobaba]

what will be formula for non uniform grid?

[@jasonbowens8369]

@@fadoobaba I believe the same. I modified the code and left that the same and I got almost exact solutions.

[@rajeebkumarmalik6388]

please send the link of pdf file of matlab code

[@naseemullahkhan9008]

i write this code and run i got singularity

[@oussamadraidi5275]

error

[@budapesteBR2012]

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

[@muhammadusama6173]

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