Unbound Shear Layer boundary condition problem

Hi, I'm very new to Matlab so bare with me, I'm trying to solve the unbounded shear layer problem for the internal section, between z=1 and -1. I cannot however enter the boundary conditon at z=1 to depend on a constant. Any help would be great.

Thanks.

[Torsten]

Please write out the problem in a mathematical notation:

Differential equation(s):

Boundary condition(s):

[Alexander Kimbley]

Hi, the differential equation is given by

y''+(a)^2 y=0, where -1<x<1.

The boundary conditions are y(-1)=0, y(1)=A. There is also an other condtion that must be satisfied which is y'(-1)=1. Noting that A is contant and not arbitrary.

Thanks.

[Torsten]

And you want to determine "a" such that all three boundary conditions are met ?

[Alexander Kimbley]

You would be correct. Obivously we can solve the problem analytically but the code will allow to explore the problem further with different conditions/criteria where one would not be able to analytically.

Answer 1

[Torsten]

function main 
  xlow = -1; xhigh = 1;
  A = 1;
  solinit = bvpinit(linspace(xlow,xhigh,4000),[1,1],1);
  sol = bvp4c(@bvp4ode, @(ya,yb,parameters)bvp4bc(ya,yb,parameters,A), solinit);
  a = sol.parameters;
  xint = linspace(xlow,xhigh,2000);
  Sxint = deval(sol,xint);
  % Analytical solution
  C1 = A/(2*sin(a));
  C2 = A/(2*cos(a));
  fun = @(x) C1*sin(a*x) + C2*cos(a*x);
  % Plot numerical vs. analytical solution
  plot(xint,Sxint(1,:),xint,fun(xint))
end 
function dydx = bvp4ode(x,y,parameters)
  dydx=[y(2); -parameters^2*y(1)];
end 
function res = bvp4bc(ya,yb,parameters,A)  
  res=[ya(1);  yb(1)-A;  ya(2)-1]; 
end

[Alexander Kimbley]

Thanks! Just one thing however, I can't seem to get the value of 'a', what code would give an output for it?

[Torsten]

But the line

a = sol.parameters 

gives the value of a.

Or what do you mean ?

[Alexander Kimbley]

Ohh I see, sorry about that! It's perfect then, thanks!

[Alexander Kimbley]

I have another problem but a lot more complex, at least for my coding skills that is. I'll attach it if you want to take a look, it'd be very appreciated. I'n the attached code I've just added the relevant ODE and BC's.

It's essentially the same problem but have introduced the c the complex eigenvalue, which is not known, B a fixed constant as well as more complex boundary conditions, I'm not sure if the last of them is needed for the problem however. R is also not known.

Thanks.

[Torsten]

As before, please state the problem in a mathematical fashion (equations and boundary conditions).

Answer 2

[Alexander Kimbley]

The ODE is given by

((x-c)^2 -B^2)(y''-ay)=2(B^2)(y-(x-c)y')/((x-c)^2), -1<x<1, B constant, a wavenumber, c complex wave speed.

Boundary conditions:

y(-1)=1,

y(1)=R,

y'(-1)=(1+c-a(1+c)^2)/(B^2 -(1+c)^2),

y'(1)=R(1-c-a(1-c)^2)/((1-c)^2 -B^2).

The last boundary condition may or may not be nessesary however.

Thanks.

[Torsten]

What constants are given ? What constants are to be determined ?

[Alexander Kimbley]

Constants given: B

Constants to be determined: c, a, R.

If the system has two many constants to be determined however, then a is probably the best to move to constants given, where a>0.

[Torsten]

function main 
  xlow = -1; xhigh = 1;
  a = 1;
  B = 0.5;
  c0 = 0.5;
  R0 = 0.5;
  y10 = 1.0;
  y20 = 0.0;
  solinit = bvpinit(linspace(xlow,xhigh,4000),[y10,y20],[c0,R0]);
  sol = bvp4c(@(x,y,p)bvp4ode(x,y,p,a,B), @(ya,yb,p)bvp4bc(ya,yb,p,a,B), solinit);
  c = sol.parameters(1)
  R = sol.parameters(2)
  xint = linspace(xlow,xhigh,2000);
  Sxint = deval(sol,xint);
  plot(xint,Sxint(1,:))
end 
function dydx = bvp4ode(x,y,p,a,B)
  c = p(1);
  R = p(2);
  dydx = [y(2); 2*B^2*(y(1)-(x-c^2)*y(2))/((x-c^2)^2*((x-c^2)^2-B^2))+a*y(1)];
end 
function res = bvp4bc(ya,yb,p,a,B)  
  c = p(1);
  R = p(2);
  res = [ya(1)-1; yb(1)-R; ya(2)-((1+c^2)-a*(1+c^2)^2)/(B^2-(1+c^2)^2); yb(2)-R*((1-c^2)-a*(1-c^2)^2)/((1-c^2)^2-B^2)]; 
end

[Alexander Kimbley]

Thanks Torsten! There just seems to be some error however, when B=0, a=0.6392, we should get that c=0 with max error being something like O(e^-5) but the output is approx 0.01 out in both complex and real parts, which leads me to think something is not quite right. I've changed the boundary conditions slightly but these do not affect the equation when B=0 either.

Is there a way to provide an inital guess for y, would this change the output?

Or could we actually make c and R variables in the function where we set them as y(3) and y(4) respectivley?

Thanks,

Alexander.

I've attached the altered code. (only the boundary conditions have changed and an error tolerance)