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

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