Loop possible? [Blasius Equation with fsolve]
显示 更早的评论
Hey ! So im pretty new to matlab, but I'm working my way into it. So now I have come across a problem, I'm not sure how to solve: I want to solve the blasius equation for a moving plate, thus different BC's than some might recognise (https://en.wikipedia.org/wiki/Blasius_boundary_layer) f''' + 1/2ff''=0, with BC's: f(0)=0, f'(0)=1 and f'(inf)=0. I know it's very easy to solve with the ODE (+shooting method) or rather the BVP command, but I need this example to apply it to a little more complex model, so I will need to do this with finite-differences(something similar to the keller box method). So since this will be boiling down to a problem formulation of a system of nonlinear equations, I figured why not use fsolve? I figured I would program a function containing all my BC's and equations, so I can just make a grid off "nodes" which will solve the equations at the corresponding node. However I seem to be stuck and or just stupid:
%%Keller box Basius with fSolve
eta0=0;
etaInf=12;
deltaEta=0.1; %stepsize
N=(etaInf-eta0)/deltaEta; %number of nodes
f=zeroes(N,1);g=ones(N,1);h=zeros(N,1); %freeing up space
anfangswert=[0,1,0.3];
%%call solver
Sol = fsolve(blasiuskeller2(f,g,h,N,deltaEta), anfangswert);
plot(eta0:edltaEta:etaInf, fval2);
---------------------------------------------------------------
%%equations mit discretiszation
function fval = blasiuskeller2(f,g,h,N,deltaEta)
% freeing up space
fval1=zeros(N-1,1);fval2=zeros(N-1,1);fval3=zeros(N-1,1);
% Define functions as F(X)=0
for i=1:N-1
fval1(i)=(f(i+1)-f(i))/deltaEta - (g(i+1)+g(i))/2;
fval2(i)=(g(i+1)-g(i))/deltaEta - (h(i+1)+h(i))/2;
fval3(i)=(h(i+1)-h(i))/deltaEta + (1/8)*(f(i+1)+f(i))*(h(i+1)+h(i));
end;
%BC's
fval(:,1)=[f(1);g(1);fval2(1)]; %initial value
for i=2:(N-1)
fval(:,i)=[fval1(i-1);fval3(i-1);fval2(i)];
end;
fval(:,N)=[fval1(N-1);1-g(N);fval2(N-1)]; %Boundary value
Pretty sure this is possible, but not sure if I'm on the right track here. Any help would be greatly appreciated!
采纳的回答
I did not check your equations, but already your call to fsolve is incorrect:
anfangswert=horzcat(zeros(1,1:N),ones(1,1:N),0.3*ones(1,1:N));
sol = fsolve(@(x)blasiuskeller2(x(1:N),x(N+1:2*N),x(2*N+1:3*N),N,deltaEta), anfangswert);
Best wishes
Torsten.
5 个评论
MaxPr
2016-8-11
Torsten
2016-8-11
Instead of
N=(etaInf-eta0)/deltaEta; %number of nodes
set
N=121;
explicitly.
Note that
number of nodes = number of intervals + 1
Best wishes
Torsten.
MaxPr
2016-8-11
Torsten
2016-8-11
anfangswert=horzcat(zeros(1,N),ones(1,N),0.3*ones(1,N));
Best wishes
Torsten.
MaxPr
2016-8-11
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Mathematics and Optimization 的更多信息
产品
Translated by 
