2D-BVP, with bvp4c loop and changing initial conditions

Hello everyone!

So I'm fairly new to matlab (meaning I did the tutorials, but aside from that I'm what you would call a noob).

I do have a fairly complicated issue: I have a PDE-BVP, which I wanna solve as an ODE-BVP in form of a loop over a bvp4c routine. So basically I have a 2D-Mesh (x and y) and I wanna solve the BVP in y direction in a loop for every x-increment.

The first and most representative part of my problem (where I'm already running into issues) has the following form: where

f'=g

g'=h

h'+f*h+2x*(h*(df/dx)-g*(dg/dx))=0

With the starting BVP's:

f(0)=0 (this should change later on to the value of f(i-1))

g(0)=1

g(inf)=0

Then I substituted the d()/dx to the finite differences formulation ()(i)-()(i-1)/deltaX, since I just want functions in the form of f(y).

The error I get right now is that the indices should be real and positiv. Which comes from my finite difference formulation for the first step. So I just tried to start the loop from i=2, however then the system has the wrong dimensions.

So I'm on the fence about my formulation.

Is there a better way to go about my problem?

Any help would be greatly appreciated!! This problem really stops me in my tracks and is haunting my for quite some time now. I already tried finite-differences together with fsolve (but that is waaaaaay to slow). I know it's solvable since I solved it with Mathematica, but for reasons I will not go into here I failed to get a useful solution later on in my model: thus Matlab ;) for better handling of the numeric issue.

Thanks in advance!!

My code sofar goes as:

function LRbvp
infinity = 10; %Eta-Inf
x0=0;
xinf=1;
deltaX=0.1;
N=((xinf-x0)/deltaX)+1; %Mesh points in x
%Initial guesses
f0=1.14267;
g0=0;
h0=0.33;

for i=1:N-1    
solinit = bvpinit(linspace(0,infinity,N),[f0 g0 h0]); %initial guess here
options = bvpset('stats','on'); 

sol = bvp4c(@(eta,f)LRode(eta,f,N,i,deltaX),@(f0,finf)LRbc(f0,finf,i)...
    ,solinit,options); %call to solver
eta = sol.x; %scalar values
f = sol.y; %colum vector of the solution
end

% --------------------------------------------------------------------------

function dfdeta = LRode(eta,f,N,i,deltaX) %Equation 37 and 39
dfdeta = [ f(2,i) %f'=g
           f(3,i) %g'=h
          -f(1,i)*f(3,i) - 2*(i/N)*(f(3,i)*((f(1,i)-f(1,i-1))/deltaX)-...
          f(2)*((f(2,i)-f(2,i-1))/deltaX))];

% --------------------------------------------------------------------------

function res = LRbc(f0,finf,i) %Boundary conditions
res = [f0(1,i) %f(0)=0
       f0(2,i) - 1 %g(0)=1
       finf(2,i)]; %g(inf)=0

[Torsten]

Maybe you could describe first in how far the model you are trying to solve differs from the usual Blasius equation.

Best wishes

Torsten.

[MaxPr]

Sure: The Blasius equation is a similar solution. With my boundary conditions I do not get a similar solution, due to the fact of a non constant heat flow (q is not constant).

So: I need to solve x & y together, I can't dismiss the x-dependency. That how I derive the form of the stream function, with f, g and h being dependent of x and y (or eta):

Answer 1

[Torsten]

bvp4c only solves the BVP on a single x-slice of your rectangular domain.

Thus you should proceed as follows:

- Solve the usual Blasius equation for x=0 and save the solution as f_old, g_old, h_old.

- Advance one step in x-direction and solve the extended Blasius equation for x=deltax. Insert fx=(f-fold)/deltax and gx=(g-gold)/deltax for the spatial derivatives in x-direction. After obtaining the solution f, g, and h for x=deltax, overwrite f_old, g_old and h_old with this new solution.

- Proceed until you reach x=xinf.

Best wishes

Torsten.

[MaxPr]

Thats exactly what I'm trying to do, but I seem to be having some issues. My question is more a "How-to" implement my thoughts into the program. I'm having issues circumventing all the errors.

So here: I'm having the issue that Matlab doesn't like my i-1 being 0 in the first step. Also when I try to do something like:

f(1,i-1)=f0;

The f0, g0 and h0 are the "old" solutions you speak of: However I seem to be unable to implement those sensibly for Matlab into the loop.

[Torsten]

The array f you get from bvp4c in LRode is a 3x1 vector (the actual solution for f,g and h in point (eta,x).

Thus you will have to supply f_old, g_old and h_old as extra (Nx1)arrays (or as a (Nx3) matrix). Something like

 sol = bvp4c(@(eta,f)LRode(eta,f,f_old,g_old,h_old,deltaX),@LRbc,solinit,options); %call to solver

should work.

The problem that arises is that you will have to find out which index in your arrays f_old, g_old and h_old corresponds to the eta value supplied by BVP4C. But you could use "interp1" to interpolate the correct f_old, g_old and h_old values.

Best wishes

Torsten.

[MaxPr]

So. I added your recommondations: I start the loop from

i=2:N

And modified my function:

function dfdeta = LRode(eta,f,N,i,deltaX,f0,g0,h0) %Equation 37 and 39
f(1,1)=f0; %f_old
g(1,1)=g0; %g_old
h(1,1)=h0; %h_old
dfdeta = [ f(2,i) %f'=g
           f(3,i) %g'=h
          -f(1,i)*f(3,i) - 2*(i/N)*(f(3,i)*((f(1,i)-f(1,i-1))/deltaX)-...
          f(2)*((f(2,i)-f(2,i-1))/deltaX))];

However now I exceed the matrix dimensions. Which makes sense because I'm calling f(2,2) in my first step. Any recommendations?

Thanks so far. I'm getting closer ^^

[Torsten]

f(2,2) is f(2) in your notation. You see that naming both "f_old" for the old "x" and "f" from BVP4C "f" will lead to inconsistencies.

Best wishes

Torsten.

[MaxPr]

Oh sure! It should be:

f(1,1)=f0; %f_old
f(2,1)=g0; %g_old
f(3,1)=h0; %h_old

Thanks for that! However the matrix dimension-error still doesn't make sense to me.

[Torsten]

I'd code it as

f_old=f0(j);
g_old=g0(j);
h_old=h0(j);

and the problem is to determine the correct "j" that corresponds to the "eta" supplied by BVP4C.

Best wishes

Torsten.