MATLAB Answers

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

20 views (last 30 days)
MaxPr
MaxPr on 13 Mar 2017
Commented: Torsten on 14 Mar 2017
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
  2 Comments
MaxPr
MaxPr on 13 Mar 2017
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):

Sign in to comment.

Accepted Answer

Torsten
Torsten on 14 Mar 2017
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.
  6 Comments
Torsten
Torsten on 14 Mar 2017
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.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!

Translated by