# fourth order Poisson iterative solver --in Matlab

I want to calculate the stream function $\psi$ starting from a velocity field $(u,v)$ (such that $u=-\frac{\partial\psi}{\partial y}$ and $v=\frac{\partial\psi}{\partial x}$). I thus calculate the vorticity $\omega$ (such that $\omega=\frac{\partial v}{\partial x}-\frac{\partial u}{\partial y}$) and then solve the Poisson equation $\nabla^2 \psi=\omega$.

To easily specify the boundary conditions, I choose a successive over-relaxation (SOR) approach.

The code seemingly runs smoothly. However, when I check if the solution is correct by differentiating $\psi$, its derivatives do not match the initial velocities $(u,v)$

Is this due to some intrinsic limitation of the SOR?

close all
clearvars
clc

%coordinates
x=1:10:600; x = repmat(x',1,50);
y=1:10:500; y = repmat(y,60,1);
dx = 10*ones(60,50); %x spacings
dy = 5*ones(60,50); %y spacings

%define velocities
u = y.^3;
v = x.^3;
figure;pcolor(x,y, u);title('u')
figure;pcolor(x,y, v);title('v')

%vorticity
figure;pcolor(x,y, dudy);title('dudy')
figure;pcolor(x,y, dvdx);title('dvdx')

omega = dvdx - dudy;

figure;pcolor(x,y, omega);title('omega');colormap(jet(10));colorbar

%iterative solver (SOR)
psi=rand(62,50); %initial guess

beta=dx./dy; %spacing ratio

%boundary conditions in y
omega(:,1:2)=0;
omega(:,end-1:end)=0;

%periodic boundary conditions in x -- append/prepend values on the other side
omega = [omega(end-1:end,:); omega; omega(1:2,:)];
beta = [beta(end-1:end,:); beta; beta(1:2,:)];

[M,N]=size(psi);

w=1.4; %relaxation factor

figure
for iter=1:1000
for ii=3:M-2
for jj=3:N-2

%%fourth order formula
psi(ii,jj) = (1-w)*psi(ii,jj) + w*...
( 16*( psi(ii+1,jj) + psi(ii-1,jj) ) - psi(ii+2,jj) - psi(ii-2,jj) +...
beta(ii,jj)^2*( -psi(ii,jj+2) + 16*(psi(ii,jj+1) + psi(ii,jj-1)) -psi(ii,jj-2) ) +...
12*dx(ii,jj)^2*omega(ii,jj) )...
/(30*(1+(beta(ii,jj)^2)));

%%second order formula
%          psi(ii,jj) = (1-w)*psi(ii,jj) + ...
%              w*( (psi(ii+1,jj)+beta(ii,jj)^2*psi(ii,jj+1)+psi(ii-1,jj)+beta(ii,jj)^2*psi(ii,jj-1)) -...
%              2*dx(ii,jj)^2*omega(ii,jj) )...
%              /(2*(1+(beta(ii,jj)^2)));

end
end

if mod(iter,10)==0
pcolor(x,y, psi(2:end-1,:))
colorbar
title(num2str(iter))
drawnow
end

end
psi = psi(2:end-1,:); %resize back

figure;pcolor(x,y, psi);title('psi')

%%%test

figure;hold on
subplot(1,2,1)
pcolor(x,y, u); colorbar; title('u orig')
subplot(1,2,2)
pcolor(x,y, -dpsidy); colorbar; title('u reconstructed')

figure;hold on
subplot(1,2,1)
pcolor(x,y, v); colorbar; title('v orig')
subplot(1,2,2)
pcolor(x,y, dpsidx); colorbar; title('v reconstructed')

• You're asking a debugging question. If you wonder whether it's due to the SOR method, have you tried to replace that part of your code by, say, a direct solver for the system of equations? – Wolfgang Bangerth Apr 4 '18 at 2:56