Finite Differences using fsolve for non linear equations

I need help writing a code that solves the boundary layer problem for fluid mechanics. I am using fsolve to solve non-linear equations (continuity and momentum equations) for the boundary layer, but I am stuck getting no solutions for my code. Any help is greatly appreciated! Here I have attempted to convert my function F (equations) and its variables, u and v, into a single vector (f and x respectively) in an attempt to solve it as I assume that fsolve takes in vector inputs. Total 2*ny*nx equations and hence unknowns to be solved, 1 for each cell.

My code:

clear all, close all, clc

nx = 10;
ny = 10;
x0 = ones(2*ny*nx,1);
k = fsolve(@blayer, x0);

function f = blayer(x)
nx = 10;
ny = 10;
L = 1;
H = 1;
dx = L/(nx-1);
dy = H/(ny-1);
P = 1;
pgrad = P/L;

v = zeros(ny,nx);
u = zeros(ny,nx);

u(:,1) = 1; % Incoming
v(:,1) = 0; % Incoming
u(1,:) = 0; % Bottom
v(1,:) = 0; % Bottom
u(ny,:) = 0; % Top
v(ny,:) = 0; % Top

F = zeros(2*ny, nx);
X = [u;v];
Xt = X';
x = Xt(:);

j = 0;
while j <= ny
    j = j+1;
    if j == 1
        for i = 1:nx
            if i == 1
                F(j,i) = (u(j,i+1)-u(j,i))/dx + (v(j+1,i)-v(j,i))/dy;
                F(j+ny,i) = u(j,i)*(u(j,i+1)-u(j,i))/dx + v(j,i)*(u(j+1,i)-u(j,i))/dy - (u(j+1,i)-2*u(j,i)+u(j,i))/dy^2 + pgrad;
            elseif i > 1 && i < nx
                F(j,i) = (u(j,i+1)-u(j,i-1))/(2*dx) + (v(j+1,i)-v(j,i))/dy;
                F(j+ny,i) = u(j,i)*(u(j,i+1)-u(j,i-1))/(2*dx) + v(j,i)*(u(j+1,i)-u(j,i))/dy - (u(j+1,i)-2*u(j,i)+u(j,i))/dy^2 + pgrad;
            elseif i == nx
                F(j,i) = (u(j,i)-u(j,i-1))/dx + (v(j+1,i)-v(j,i))/dy;
                F(j+ny,i) = u(j,i)*(u(j,i)-u(j,i-1))/dx + v(j,i)*(u(j+1,i)-u(j,i))/dy - (u(j+1,i)-2*u(j,i)+u(j,i))/dy^2 + pgrad;
            end
        end
    end
    if j > 1 && j < ny
        for i = 1:nx
            if i == 1
                F(j,i) = (u(j,i+1)-u(j,i))/dx + (v(j+1,i)-v(j-1,i))/(2*dy);
                F(j+ny,i) = u(j,i)*(u(j,i+1)-u(j,i))/dx + v(j,i)*(u(j+1,i)-u(j-1,i))/(2*dy) - (u(j+1,i)-2*u(j,i)+u(j-1,i))/dy^2 + pgrad;
            elseif i > 1 && i < nx
                F(j,i) = (u(j,i+1)-u(j,i-1))/(2*dx) + (v(j+1,i)-v(j-1,i))/(2*dy);
                F(j+ny,i) = u(j,i)*(u(j,i+1)-u(j,i-1))/(2*dx) + v(j,i)*(u(j+1,i)-u(j-1,i))/(2*dy) - (u(j+1,i)-2*u(j,i)+u(j-1,i))/dy^2 + pgrad;
            elseif i == nx
                F(j,i) = (u(j,i)-u(j,i-1))/dx + (v(j+1,i)-v(j-1,i))/(2*dy);
                F(j+ny,i) = u(j,i)*(u(j,i)-u(j,i-1))/dx + v(j,i)*(u(j+1,i)-u(j-1,i))/(2*dy) - (u(j+1,i)-2*u(j,i)+u(j-1,i))/dy^2 + pgrad;
            end
        end
    end
    if j == ny
        for i = 1:nx
            if i == 1
                F(j,i) = (u(j,i+1)-u(j,i))/dx + (v(j,i)-v(j-1,i))/dy;
                F(j+ny,i) = u(j,i)*(u(j,i+1)-u(j,i))/dx + v(j,i)*(u(j,i)-u(j-1,i))/dy - (u(j,i)-2*u(j,i)+u(j-1,i))/dy^2 + pgrad;
            elseif i > 1 && i < nx
                F(j,i) = (u(j,i+1)-u(j,i-1))/(2*dx) + (v(j,i)-v(j-1,i))/dy;
                F(j+ny,i) = u(j,i)*(u(j,i+1)-u(j,i-1))/(2*dx) + v(j,i)*(u(j,i)-u(j-1,i))/dy - (u(j,i)-2*u(j,i)+u(j-1,i))/dy^2 + pgrad;
            elseif i == nx
                F(j,i) = (u(j,i)-u(j,i-1))/dx + (v(j,i)-v(j-1,i))/dy;
                F(j+ny,i) = u(j,i)*(u(j,i)-u(j,i-1))/dx + v(j,i)*(u(j,i)-u(j-1,i))/dy - (u(j,i)-2*u(j,i)+u(j-1,i))/dy^2 + pgrad;
            end
        end
    end
end

Ft = F';
f = Ft(:);
end

[Ellson Chow]

I have changed my initial guess to u = 1 and v = 0;

x0 = zeros(2*ny*nx,1);

x0(1:ny*nx) = 1;

but still does not work.

Answer 1

[Alan Stevens]

Do you need fsolve at all? Doesn't your blayer function do the solving itself? Couldn't you simply have

nx = 10;
ny = 10;
L = 1;
H = 1;
dx = L/(nx-1);
dy = H/(ny-1);
P = 1;
pgrad = P/L;
v = zeros(ny,nx);
u = zeros(ny,nx);
u(:,1) = 1; % Incoming
v(:,1) = 0; % Incoming
u(1,:) = 0; % Bottom
v(1,:) = 0; % Bottom
u(ny,:) = 0; % Top
v(ny,:) = 0; % Top
F = zeros(2*ny, nx);
% X = [u;v];
% Xt = X';
% x = Xt(:);
j = 0;
while j <= ny
    j = j+1;
    if j == 1
        for i = 1:nx
            if i == 1
                F(j,i) = (u(j,i+1)-u(j,i))/dx + (v(j+1,i)-v(j,i))/dy;
                F(j+ny,i) = u(j,i)*(u(j,i+1)-u(j,i))/dx + v(j,i)*(u(j+1,i)-u(j,i))/dy - (u(j+1,i)-2*u(j,i)+u(j,i))/dy^2 + pgrad;
            elseif i > 1 && i < nx
                F(j,i) = (u(j,i+1)-u(j,i-1))/(2*dx) + (v(j+1,i)-v(j,i))/dy;
                F(j+ny,i) = u(j,i)*(u(j,i+1)-u(j,i-1))/(2*dx) + v(j,i)*(u(j+1,i)-u(j,i))/dy - (u(j+1,i)-2*u(j,i)+u(j,i))/dy^2 + pgrad;
            elseif i == nx
                F(j,i) = (u(j,i)-u(j,i-1))/dx + (v(j+1,i)-v(j,i))/dy;
                F(j+ny,i) = u(j,i)*(u(j,i)-u(j,i-1))/dx + v(j,i)*(u(j+1,i)-u(j,i))/dy - (u(j+1,i)-2*u(j,i)+u(j,i))/dy^2 + pgrad;
            end
        end
    end
    if j > 1 && j < ny
        for i = 1:nx
            if i == 1
                F(j,i) = (u(j,i+1)-u(j,i))/dx + (v(j+1,i)-v(j-1,i))/(2*dy);
                F(j+ny,i) = u(j,i)*(u(j,i+1)-u(j,i))/dx + v(j,i)*(u(j+1,i)-u(j-1,i))/(2*dy) - (u(j+1,i)-2*u(j,i)+u(j-1,i))/dy^2 + pgrad;
            elseif i > 1 && i < nx
                F(j,i) = (u(j,i+1)-u(j,i-1))/(2*dx) + (v(j+1,i)-v(j-1,i))/(2*dy);
                F(j+ny,i) = u(j,i)*(u(j,i+1)-u(j,i-1))/(2*dx) + v(j,i)*(u(j+1,i)-u(j-1,i))/(2*dy) - (u(j+1,i)-2*u(j,i)+u(j-1,i))/dy^2 + pgrad;
            elseif i == nx
                F(j,i) = (u(j,i)-u(j,i-1))/dx + (v(j+1,i)-v(j-1,i))/(2*dy);
                F(j+ny,i) = u(j,i)*(u(j,i)-u(j,i-1))/dx + v(j,i)*(u(j+1,i)-u(j-1,i))/(2*dy) - (u(j+1,i)-2*u(j,i)+u(j-1,i))/dy^2 + pgrad;
            end
        end
    end
    if j == ny
        for i = 1:nx
            if i == 1
                F(j,i) = (u(j,i+1)-u(j,i))/dx + (v(j,i)-v(j-1,i))/dy;
                F(j+ny,i) = u(j,i)*(u(j,i+1)-u(j,i))/dx + v(j,i)*(u(j,i)-u(j-1,i))/dy - (u(j,i)-2*u(j,i)+u(j-1,i))/dy^2 + pgrad;
            elseif i > 1 && i < nx
                F(j,i) = (u(j,i+1)-u(j,i-1))/(2*dx) + (v(j,i)-v(j-1,i))/dy;
                F(j+ny,i) = u(j,i)*(u(j,i+1)-u(j,i-1))/(2*dx) + v(j,i)*(u(j,i)-u(j-1,i))/dy - (u(j,i)-2*u(j,i)+u(j-1,i))/dy^2 + pgrad;
            elseif i == nx
                F(j,i) = (u(j,i)-u(j,i-1))/dx + (v(j,i)-v(j-1,i))/dy;
                F(j+ny,i) = u(j,i)*(u(j,i)-u(j,i-1))/dx + v(j,i)*(u(j,i)-u(j-1,i))/dy - (u(j,i)-2*u(j,i)+u(j-1,i))/dy^2 + pgrad;
            end
        end
    end
end

surf(F)

or am I missing something?

[Ellson Chow]

Hi thanks for answering. I am attempting to solve for u(i,j) and v(i,j) in a grid of nx*ny cells. As the momentum equation has u(i,j)*(du/dx) and v(i,j)*(du/dy), this becomes a non-linear equation when I expand du/dx and du/dy using central difference method. Thus, I assume that I need to use fsolve due to its non-linearity. I tried your code and it does show me the values for u(i,j) and v(i,j) in a nx*ny grid with its initial values, but not solved (basically everything else is still 0).

[Alan Stevens]

So you want to find the values of u and v that make F zero everywhere? In that case you need to pass u and v to blayer, not x (the values of which are fixed). i.e. you want fsolve to find values of u and v that set all the Fs to zero.

[Abdallah Daher]

Were you able it out? I am currently trying to write an explicit finite difference scheme for a flat plate boundary layer