Upwind, Central Differencing, and Upwind 2nd Order Solution to PDE

12 views (last 30 days)
Hello,
I am trying to find a solution to this PDE:
u*+v=
u=sin(45)
v=cos(45)
*10^-6
The solution space is a 1X1 square with the following boundary conditions:
@
@
@ x
@ x=1
@ y=1
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

Answers (0)

Categories

Find more on Labeling, Segmentation, and Detection in Help Center and File Exchange

Products


Release

R2018b

Community Treasure Hunt

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

Start Hunting!