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

12 次查看（过去 30 天）
Daria Nikitaeva
Daria Nikitaeva 2019-10-10
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

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Labeling, Segmentation, and Detection 的更多信息

标签

产品


版本

R2018b

Community Treasure Hunt

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

Start Hunting!

Translated by