How can I solve BVP with multiple boundary conditions?

Question

I am facing some problems in solving my differential equations using the forward difference method. I have two regions (0 to t and t to 1, and t is any value from 0 to 1). In first region (0 to t), I have 4 differential equations and in the second region (t to 1) I have 2 differential equations, total I have six boundary conditions (BC), 2 BC in the first region, 4 BC at an interface (t) and 1 BC in the second region.

I followed the documentation in MathWorks to solve BVP with multiple boundary conditions (https://www.mathworks.com/help/matlab/math/solve-bvp-with-multiple-boundary-conditions.html)

When I run my code I am getting an error message, not enough input arguments. I think the problem lies in the region. BVP solver is not getting the required region to solve the equations. Please, anyone, guide me on how to overcome this problem? it will be appreciated.

t=0.9;    
x = [0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 t t 0.95 1];
solinit = bvpinit(x,[zeros(1,90)]);
sol = bvp5c(@ex1ode,@ex1bc,solinit);
y = deval(sol,x);
plot(sol.x,sol.y(1,:));
function dydx = ex1ode(x,y,region)
h = 0.0001; t = 0.9; k = 100; Bi = 0.01; Da = 0.0001;
A0 = 100; A1 = 0.944545454545455; A2 = -0.4445454545454545; A3 = 7.3001203377370034e-43; A6 = 0.0101;
E1 = 0.002414161912718531; E2 = 3.55998187608043e-45 ; E3 = 0.00476100894331544; E4 = 0; E5 = -0.816974047611677;
D1 = 6.785528299098551; D2 = -10.83936751884141; D3 = 7.676961726023331; D4 = -2.03191960987143; D5 = -16.67685194071544;
uc=((-(x^2)/2)+A1*x+A2)/(-0.166667+(0.5*A1)+A2-(A2*t)+(Da*t)-(0.5*A1*(t^2))+0.166667*(t^3)+ 
(A3*sinh(A0*t)/A0));
up =(A0*(A3*cosh(A0*x)+Da))/(A0(-0.166667+A2-(A2*t)+(Da*t)+0.166667*(t^3)+A1*(0.5-0.5* 
(t^2)))+A3*sinh(A0*t));
yini=-E1*(x^2)-E2*cosh(A0*x)-E3*cosh(sqrt(A6)*x)-E4*cosh(2*A0*x)-E5;
tini=-D1*x-D2*(x^2)-D3*(x^3)-D4*(x^4)-D5;
n=90;
dydx=zeros(6,1);
for i=1:6:n
    switch region
        case 1  % x in [0 t]
            dydx(i)=y(i+1);
            dydx(i+1)=up*((y(i)-yini)/h)-(Bi./k)*(y(i+2)-y(i));
            dydx(i+2)=y(i+3);
            dydx(i+3)=Bi*(y(i+2)-y(i));
            yini=y(i);
        case 2  % x in [t 1]
            dydx(i+4)=y(i+5);
            dydx(i+5)=uc*((y(i+4)-tini)/h);
            tini=y(i+4);
end
end
end

function res=ex1bc(YL,YR)
n=90;k =0.01;
res=zeros(n,1);
for i=1:6:n
    res(i)=YL(i+1,1);
    res(i+1)=YL(i+3,1);
    res(i+2)=YR(i,1)-YR(i+2,1);
    res(i+3)=YR(i+2,1)-YL(i+4,2);
    res(i+4)=k*YL(i+5,2)-k*YR(i+1,1)-YR(i+3,1);
    res(i+5)=YR(i+5,2)-1;
end
end

"I am having some bug in my code" Please provide more information about what you are experiencing. Do you get wrong results? An error message? Does your computer hang? Does it erase your hard drive? Please read [ask], then [edit] your code accordingly. — Cris Luengo, Apr 03 '20 at 19:43
Does your system have dimension 6=4+2 or 90=15*(4+2)? Where is the region determined? Why do you not use the method provided by Matlab for 3 or more "boundary" points? Can you produce a simplified test case for a MWE? — Lutz Lehmann, Apr 04 '20 at 06:35
I have used the same method as in Matlab, they defined the region the same as this. (https://www.mathworks.com/help/matlab/math/solve-bvp-with-multiple-boundary-conditions.html). I have dimension 90=15*(4+2). What does MWE mean? — MIRZA FARRUKH BAIG, Apr 04 '20 at 06:59
It means a minimal working example, reducing dimensions and simplifying equation (thus changing them) while keeping the principal structure and so that the error is still produced. — Lutz Lehmann, Apr 04 '20 at 08:46
You are missing the continuity boundary conditions, the return array should have dimension `2*n`. You should initialize the derivatives vector with the actual dimension or you miss the trailing zeros in region 1. — Lutz Lehmann, Apr 04 '20 at 09:00
Do you mean continuity BC at the interface? If yes then I have Bc at the interface representing in res(i+2), res(i+3) and res(i+4). — MIRZA FARRUKH BAIG, Apr 04 '20 at 09:46
I mean where `YR(k,1)==YL(k,2)` for all `k=1:n`. This ensures continuity of the solution components. You have essentially 2 boundary value problems, and each needs `n` boundary conditions. — Lutz Lehmann, Apr 04 '20 at 10:09
This means I need to add more boundary conditions to ensure continuity at the interface. — MIRZA FARRUKH BAIG, Apr 04 '20 at 11:17
@LutzLehmann can you edit my code? So I can understand easily. — MIRZA FARRUKH BAIG, Apr 06 '20 at 10:36
I do not have matlab and octave has no usable bvp4c implementation. You would get a much easier to understand mockup problem if you used something like oscillators or dampened spring equations in both regions. The matlab example is somewhat unfortunate as they put the trivial (but necessary) continuity conditions in the middle of the boundary conditions. — Lutz Lehmann, Apr 06 '20 at 10:43
thanks for your suggestions and help. Ok, I will try to solve but still have many confusions at the interface boundary. — MIRZA FARRUKH BAIG, Apr 06 '20 at 10:53

How can I solve BVP with multiple boundary conditions?

0 Answers0