Upwind, Central Differencing, and Upwind 2nd Order Solution to PDE
5 views (last 30 days)
Show older comments
Hello,
I am trying to find a solution to this PDE:
u*
+v
=![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/242102/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/242100/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/242101/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/242102/image.png)
u=sin(45)
v=cos(45)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/242103/image.png)
The solution space is a 1X1 square with the following boundary conditions:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/242104/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/242105/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/242106/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/242107/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/242108/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/242109/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/242110/image.png)
The goal is to compare central differencing, upwind, and upwind 2nd order solutions for ϕ at y=0.5 at the following grid sizes: 11X11, 21X21, and 41X41.
I wrote the following code, however, my professor says that it's incorrect. I tried putting it into the PDE Modeler in MATLAB and it gives me the same result as my upwind solution.
I'm lost at this point. Can someone please help?
Thank you!
clc
clear all
close all
% compare CD2 UD1, and UD2 on the same figure
% each figure will be for one resolution
% Need to use a mirror image for the north and east sides (the derivatives
% are equal to zero)
alpha=10^-6;
theta=45;
u=cosd(theta);
v=sind(theta);
grid_size=[11,21,41];
%grid_size=[11];
for k=1:length(grid_size)
grid=linspace(0,1,grid_size(k));
resolution(k)=grid(2)-grid(1);
dx=resolution(k);
dy=resolution(k);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%CD2%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A=zeros(grid_size(k)^2);
%create general solution pattern
%main diagonal
for i=1:grid_size(k)^2
A(i,i)=2*(alpha/dx^2+alpha/dy^2);
B(i)=0;
end
%horizontals
for i=1:grid_size(k)-1
for j=1:grid_size(k)
A(i+(j-1)*grid_size(k),i+(j-1)*grid_size(k)+1)=u/(2*dx)-alpha/dx^2;
A(i+(j-1)*grid_size(k)+1,i+(j-1)*grid_size(k))=-u/(2*dx)-alpha/dx^2;
end
end
%verticals
for i=1:grid_size(k)
for j=1:grid_size(k)-1
A(i+(j-1)*grid_size(k),i+j*grid_size(k))=v/(2*dx)-alpha/dy^2;
A(i+j*grid_size(k),i+(j-1)*grid_size(k))=-v/(2*dx)-alpha/dy^2;
end
end
%establish boundary conditions
%west boundary
i=1;
for j=1:grid_size(k):grid_size(k)^2
A(j,:)=0;
A(j,j)=1;
if grid(i)<=0.2
B(j)=1;
end
i=i+1;
end
%south boundary
for i=1:grid_size(k)
A(i,:)=0;
A(i,i)=1;
B(i)=1;
end
%north boundary
%main diagonal
for i=grid_size(k)^2-grid_size(k)+1:grid_size(k)^2
A(i,:)=0;
A(i,i)=2*alpha/dx^2;
B(i)=0;
end
%horizontals
for i=1:grid_size(k)-1
for j=grid_size(k)
A(i+(j-1)*grid_size(k),i+(j-1)*grid_size(k)+1)=u/(2*dx)-alpha/dx^2;
A(i+(j-1)*grid_size(k)+1,i+(j-1)*grid_size(k))=-u/(2*dx)-alpha/dx^2;
end
end
%east boundary
%main diagonal
for i=grid_size(k):grid_size(k):grid_size(k)^2
A(i,:)=0;
A(i,i)=2*alpha/dx^2;
B(i)=0;
end
%verticals
for i=grid_size(k)
for j=1:grid_size(k)-1
A(i+(j-1)*grid_size(k),i+j*grid_size(k))=v/(2*dx)-alpha/dy^2;
A(i+j*grid_size(k),i+(j-1)*grid_size(k))=-v/(2*dx)-alpha/dy^2;
end
end
%A=sparse(A);
%B=sparse(B);
%dA=decomposition(A);
%phi_vector=dA\B';
%phi_vector=A\B';
%phi_vector=mldivide(A,B');
phi_vector=linsolve(A,B');
%phi_vector=inv(A)*B';
for i=1:length(grid)
for j=1:length(grid)
phi_CD2(i,j)=phi_vector(length(grid)*(i-1)+j);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%UD1%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear A B phi_vector
A=zeros(grid_size(k)^2);
%create general solution pattern
%main diagonal
for i=1:grid_size(k)^2
A(i,i)=u/dx+v/dy+2*(alpha/dx^2+alpha/dy^2);
B(i)=0;
end
%horizontals
for i=1:grid_size(k)-1
for j=1:grid_size(k)
A(i+(j-1)*grid_size(k),i+(j-1)*grid_size(k)+1)=-alpha/dx^2;
A(i+(j-1)*grid_size(k)+1,i+(j-1)*grid_size(k))=-u/(2*dx)-alpha/dx^2;
end
end
%verticals
for i=1:grid_size(k)
for j=1:grid_size(k)-1
A(i+(j-1)*grid_size(k),i+j*grid_size(k))=-alpha/dy^2;
A(i+j*grid_size(k),i+(j-1)*grid_size(k))=-v/(2*dx)-alpha/dy^2;
end
end
%establish boundary conditions
%west boundary
i=1;
for j=1:grid_size(k):grid_size(k)^2
A(j,:)=0;
A(j,j)=1;
if grid(i)<=0.2
B(j)=1;
end
i=i+1;
end
%south boundary
for i=1:grid_size(k)
A(i,:)=0;
A(i,i)=1;
B(i)=1;
end
%north boundary
%main diagonal
for i=grid_size(k)^2-grid_size(k)+1:grid_size(k)^2
A(i,:)=0;
A(i,i)=u/dx+2*alpha/dx^2;
B(i)=0;
end
%horizontals
for i=1:grid_size(k)-1
for j=grid_size(k)
A(i+(j-1)*grid_size(k),i+(j-1)*grid_size(k)+1)=-alpha/dx^2;
A(i+(j-1)*grid_size(k)+1,i+(j-1)*grid_size(k))=-u/(2*dx)-alpha/dx^2;
end
end
%east boundary
%main diagonal
for i=grid_size(k):grid_size(k):grid_size(k)^2
A(i,:)=0;
A(i,i)=v/dy+2*alpha/dx^2;
B(i)=0;
end
%verticals
for i=grid_size(k)
for j=1:grid_size(k)-1
A(i+(j-1)*grid_size(k),i+j*grid_size(k))=-alpha/dy^2;
A(i+j*grid_size(k),i+(j-1)*grid_size(k))=-v/(2*dx)-alpha/dy^2;
end
end
%A=sparse(A);
%B=sparse(B);
%dA=decomposition(A);
%phi_vector=dA\B';
%phi_vector=mldivide(A,B');
%phi_vector=A\B';
%phi_vector=inv(A)*B';
phi_vector=linsolve(A,B');
for i=1:length(grid)
for j=1:length(grid)
phi_UD1(i,j)=phi_vector(length(grid)*(i-1)+j);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%UD2%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear A B phi_vector
A=zeros(grid_size(k)^2);
%create general solution pattern
%main diagonal
for i=1:grid_size(k)^2
A(i,i)=3*u/(2*dx)+3*v/(2*dy)+2*(alpha/dx^2+alpha/dy^2);
B(i)=0;
end
%horizontals
for i=1:grid_size(k)-2
for j=1:grid_size(k)
A(i+(j-1)*grid_size(k),i+(j-1)*grid_size(k)+1)=-alpha/dx^2;
A(i+(j-1)*grid_size(k)+1,i+(j-1)*grid_size(k))=-3*u/(2*dx)-alpha/dx^2;
A(i+(j-1)*grid_size(k)+2,i+(j-1)*grid_size(k))=3*u/(2*dx);
end
end
%verticals
for i=1:grid_size(k)
for j=2:grid_size(k)-1
A(i+(j-1)*grid_size(k),i+j*grid_size(k))=-alpha/dy^2;
A(i+j*grid_size(k),i+(j-1)*grid_size(k))=-3*v/(2*dx)-alpha/dy^2;
A(i+j*grid_size(k),i+(j-2)*grid_size(k))=3*v/(2*dx);
end
end
%establish boundary conditions
%west boundary
i=1;
for j=1:grid_size(k):grid_size(k)^2
A(j,:)=0;
A(j,j)=1;
if grid(i)<=0.2
B(j)=1;
end
i=i+1;
end
i=1;
for j=2:grid_size(k):grid_size(k)^2
A(j,:)=0;
A(j,j)=1;
if grid(i)<=0.2
B(j)=1;
end
i=i+1;
end
%south boundary
for i=1:grid_size(k)*2
A(i,:)=0;
A(i,i)=1;
B(i)=1;
end
%north boundary
%main diagonal
for i=grid_size(k)^2-grid_size(k)+1:grid_size(k)^2
A(i,:)=0;
A(i,i)=3*u/(2*dx)+2*alpha/dx^2;
B(i)=0;
end
%horizontals
for i=1:grid_size(k)-2
for j=grid_size(k)
A(i+(j-1)*grid_size(k),i+(j-1)*grid_size(k)+1)=-alpha/dx^2;
A(i+(j-1)*grid_size(k)+1,i+(j-1)*grid_size(k))=-3*u/(2*dx)-alpha/dx^2;
A(i+(j-1)*grid_size(k)+2,i+(j-1)*grid_size(k))=3*u/(2*dx);
end
end
%east boundary
%main diagonal
for i=grid_size(k):grid_size(k):grid_size(k)^2
A(i,:)=0;
A(i,i)=3*v/(2*dy)+2*alpha/dx^2;
B(i)=0;
end
%verticals
for i=grid_size(k)
for j=2:grid_size(k)-1
A(i+(j-1)*grid_size(k),i+j*grid_size(k))=-alpha/dy^2;
A(i+j*grid_size(k),i+(j-1)*grid_size(k))=-3*v/(2*dx)-alpha/dy^2;
A(i+j*grid_size(k),i+(j-2)*grid_size(k))=3*v/(2*dx);
end
end
%A=sparse(A);
%B=sparse(B);
%dA=decomposition(A);
%phi_vector=dA\B';
%phi_vector=mldivide(A,B');
%phi_vector=A\B';
%phi_vector=inv(A)*B';
phi_vector=linsolve(A,B');
for i=1:length(grid)
for j=1:length(grid)
phi_UD2(i,j)=phi_vector(length(grid)*(i-1)+j);
end
end
figure(k)
subplot(2,2,1)
surf(grid,grid,phi_CD2)
shading interp
xlabel('x')
ylabel('y')
title(['CD2',' ',num2str(length(grid)),'X',num2str(length(grid))])
zlim([0 1])
subplot(2,2,2)
surf(grid,grid,phi_UD1)
shading interp
xlabel('x')
ylabel('y')
title(['UD1',' ',num2str(length(grid)),'X',num2str(length(grid))])
zlim([0 1])
subplot(2,2,3)
surf(grid,grid,phi_UD2)
shading interp
xlabel('x')
ylabel('y')
title(['UD2',' ',num2str(length(grid)),'X',num2str(length(grid))])
zlim([0 1])
subplot(2,2,4)
plot(grid,phi_CD2((grid_size(k)-1)/2,:))
hold on
plot(grid,phi_UD1((grid_size(k)-1)/2,:))
plot(grid,phi_UD2((grid_size(k)-1)/2,:))
title(['\phi at y=0.5',' ',num2str(length(grid)),'X',num2str(length(grid))])
legend('CD2','UD1','UD2','Location','best')
xlabel('x')
ylabel('\phi')
end
0 Comments
Answers (0)
See Also
Categories
Find more on General PDEs in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!