How to solve equation by iteration

Hello,

Im trying to solve this equaiton below (for Psi(i,j)) but im not sure what im missing. I should be getting a 220x40 matrix with various values for Psi but im not sure how to go about doing this.

clear
clc
Vinf = 1;
c = 1;
L = 11*c;
H = 4*c;
N = 220;
deltax = L/N;
M = 40;
deltay = H/M;
epsilon = 10^(-6);
 
yi=linspace(0,L,M);
xi=linspace(0,H,N);
z=zeros(M,N);
for i = 1:N-1
    for j = 1:M-1
        x(i,j) = i*deltax;
        y(i,j) = j*deltay;
        Psi(i,j) = ((deltay^2)/(2*((deltax^2)+(deltay^2))))*(Psi(i+1,j)+Psi(i-1,j))+((deltax^2)/(2*((deltax^2)+(deltay^2))))*(Psi(i,j+1)+Psi(i,j-1));
    end
end
%boundary conditions
Psi(0,j) = Vinf*deltay*j;
Psi(i,M+1) = Vinf*H;
Psi(N,j) = Psi(N-1,j);
Psi(i,1) = 0;

[Torsten]

Your equations can't be solved directly for Psi because Psi appears on both sides of the equation.

But your equations define a linear system for Psi.

Determine A such that your system is of the form

A*Psi = b

and solve for Psi as

Psi = A\b.

Or if you are told to use Jacobi iteration, apply a scheme like

Psi_new(i,j) = function(Psi_old(i,j),Psi_old(i+1,j),Psi_old(i-1,j),Psi_old(i,j+1),Psi_old(i,j+1))

Thus you will not only need one matrix Psi, but two matrices Psi_old and Psi_new.

Answer 1

[Star Strider]

One problem is that ‘Psi’ is used before it is defined. Preallocating fixes that, however it may not solve the larger problems. There were also indexing problems that were fixed by beginning both loops at 2 rather than 1. I added the surf plot to visualise the result.

Since I have no idea what the objective is or what the mathematical relations are, I have no other suggestions.

Vinf = 1;
c = 1;
L = 11*c;
H = 4*c;
N = 220;
deltax = L/N;
M = 40;
deltay = H/M;
epsilon = 10^(-6);
 
yi=linspace(0,L,M);
xi=linspace(0,H,N);
z=zeros(M,N);
Psi = zeros(N,M);                                       % Preallocate
for i = 2:N-1
    for j = 2:M-1
        x(i,j) = i*deltax;
        y(i,j) = j*deltay;
        Psi(i,j) = ((deltay^2)/(2*((deltax^2)+(deltay^2))))*(Psi(i+1,j)+Psi(i-1,j))+((deltax^2)/(2*((deltax^2)+(deltay^2))))*(Psi(i,j+1)+Psi(i,j-1));
    end
end
%boundary conditions
Psi(1,j) = Vinf*deltay*j;
Psi(i,M+1) = Vinf*H;
Psi(N,j) = Psi(N-1,j);
Psi(i,1) = 0;

figure
surf(Psi, 'EdgeColor','none')
grid on

.