bvp4c (Compressible boundary layer)

Hi, i hope someone can help me. I want to find the solution to the compressible boundary layer equations, this problem is part of my thesis project, but I'm running into some problems.

To make it easier to read, here is the coupled set of ODE that I'm trying to solve;

y(4)=1/C*(-C'*y(3)-y(1)*y(3))

y(7)=1/C*(-P*y(1)*y(6)-P*C*(g-1)M^2*y(3)^2-C'*y(6))

Where C and C' are;

C=(y(5)^0.5*(h+c*s)/(y(5)*h+c*s))

C'=(y(6)*((T/s)+1)*(1-(T/s)*y(5))/(2*y(5)^0.5*(y(5)*(T/s)+1)^2))

My code

function Compressible_Boundary_Layer
h=301350; c=1004.5; s=120; T=300; P=1; g=1.4; M=2.9;
solinit=bvpinit(linspace(0,8,9),@Compressibleguess);
sol=bvp4c(@Compressible,@Compressiblebc,solinit,[],h,c,s,T,P,M);
x=linspace(0,8,100);
y=deval(sol,x);
plot(x,y(2,:))
function dydx=Compressible(x,y,h,c,s,T,P,M)
dydx=[y(2);
y(3);
1/(y(5)^0.5*(h+c*s)/(y(5)*h+c*s))*(-(y(6)*((T/s)+1)*(1-(T/s)*y(5))/(2*y(5)^0.5*(y(5)*(T/s)+1)^2))*y(3)-y(1)*y(3));
y(6);
1/(y(5)^0.5*(h+c*s)/(y(5)*h+c*s))*(-P*y(1)*y(6)-P*(y(5)^0.5*(h+c*s)/(y(5)*h+c*s))*(g-1)*M^2*y(3)^2-(y(6)*((T/s)+1)*(1-(T/s)*y(5))/(2*y(5)^0.5*(y(5)*(T/s)+1)^2))*y(6))];
function res=Compressiblebc(ya,yb,h,c,s,T,P,M)
res=[ya(1);ya(2);yb(2)-1;yb(5)-1;ya(6)];
function y=Compressibleguess(x)
y(1)=x;
y(2)=x^0.5;
y(3)=5-x;
y(5)=-x^2+2;
y(6)=-x;

When i run this i get two errors;

??? Error using ==> bvp4c The boundary condition function should return a column vector of length 6

Error in ==> Compressible_Boundary_Layer at 4 sol=bvp4c(@Compressible,@Compressiblebc,solinit,[],h,c,s,T,P,M);

Hope someone can help, thanks in advance.

Ayden

Answer 1

[Walter Roberson]

Before assigning to y(1), insert

y = zeros(6,1);

By default arrays accessed with a single index are created as row vectors rather than column vectors.

[Ayden]

Thanks for the suggestion. Unfortunately I still end up with the same error when i insert

y=zeros(6,1);

It appears to be asking me for 6 boundary conditions when i only have 5 first order equations.

[Walter Roberson]

I just noticed that you do not assign anything to y(4) in your guess function?

[Walter Roberson]

Might I suggest that you only want your guess function to create y(1) to y(5) ?

[Ayden]

I though I didn't need to specify the intial guess function for y(4) because it depends on the guess of y(1) and y(3) (which i have already specified), I also didn't specify y(7) because of this.

y(4)=1/C*(-C'*y(3)-y(1)*y(3))

y(7)=1/C*(-P*y(1)*y(6)-P*C*(g-1)M^2*y(3)^2-C'*y(6))

Do I still need to assign guess functions to y(4) and y(7), cause i already have 5 first order DE, with 5 boundary conditions and 5 initial guess functions at the moment.

Answer 2

[Ayden]

Thanks for your help so far Walter, i managed to fix my initial problem by placing the guess functions into a matrix. I dont encounter that error anymore.

I also corrected some errors in my inital code relating to the unknown parameters specified in bvp4c. I removed them since they are just known parameters for the DE functions.

my new code

function Compressible_Boundary_Layer
solinit=bvpinit(linspace(0,8,9),@guess);
sol=bvp4c(@Compressible,@Compressiblebc,solinit);
x=linspace(0,8,100);
y=deval(sol,x);
plot(x,y(2,:))
%---------------------------------------------------
function dydx=Compressible(x,y)
h=301350; c=1004.5; s=120; T=300; P=1; g=1.4; M=2.9;
%y = zeros(6,1);
dydx=[y(2);
    y(3);
    1/(y(5)^0.5*(h+c*s)/(y(5)*h+c*s))*(-1*(y(6)*((T/s)+1)*(1-(T/s)*y(5))/(2*y(5)^0.5*(y(5)*(T/s)+1)^2))*y(3)-y(1)*y(3));
    y(6);
    1/(y(5)^0.5*(h+c*s)/(y(5)*h+c*s))*(-P*y(1)*y(6)-P*(y(5)^0.5*(h+c*s)/(y(5)*h+c*s))*(g-1)*M^2*y(3)^2-(y(6)*((T/s)+1)*(1-(T/s)*y(5))/(2*y(5)^0.5*(y(5)*(T/s)+1)^2))*y(6))];
%---------------------------------------------------
function res=Compressiblebc(ya,yb)
res=[ya(1);ya(2);yb(2)-1;yb(5)-1;ya(6)];
%---------------------------------------------------
function y=guess(x)
y=[x;x^0.5;5-x;exp(-x);-4/5/exp(2*(x-2)^2)];

This time round, it appears that i have some issues with my boundary conditions, this is the error that I'm getting;

??? Attempted to access ya(6); index out of bounds because numel(ya)=5.

Any ideas on what is causing this. Thanks in advance.