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

10 次查看（过去 30 天）
MaxPr
MaxPr 2017-3-13
评论： Torsten 2017-3-14
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 个评论
Torsten
Torsten 2017-3-13
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
MaxPr 2017-3-13
编辑：MaxPr 2017-3-13
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):

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2017-3-14
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 个评论
显示 4更早的评论
MaxPr
MaxPr 2017-3-14
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
Torsten 2017-3-14
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.

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Boundary Value Problems 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by