Unbound Shear Layer boundary condition problem

3 次查看（过去 30 天）
Alexander Kimbley
编辑： Alexander Kimbley 2019-4-18
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.
  4 个评论
显示 2更早的评论
Torsten
Torsten 2019-4-10
And you want to determine "a" such that all three boundary conditions are met ?
Alexander Kimbley
Alexander Kimbley 2019-4-10
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.

请先登录，再进行评论。

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

采纳的回答

Torsten
Torsten 2019-4-10
编辑：Torsten 2019-4-10
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
  5 个评论
显示 3更早的评论
Alexander Kimbley
Alexander Kimbley 2019-4-12
编辑：Alexander Kimbley 2019-4-12
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
Torsten 2019-4-15
As before, please state the problem in a mathematical fashion (equations and boundary conditions).

请先登录，再进行评论。

更多回答（1 个）

Alexander Kimbley
Alexander Kimbley 2019-4-16
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.
  4 个评论
显示 2更早的评论
Torsten
Torsten 2019-4-16
编辑：Torsten 2019-4-16
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
Alexander Kimbley 2019-4-18
编辑：Alexander Kimbley 2019-4-18
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)

请先登录，再进行评论。

类别

在 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