BVP5C multipoint boundary conditions sequence

Hello everyone,

I'm trying to solve a system of two coupled second order matrix (2x2) differential equations with multiple boundaries; 4 to be exact. While my code is working, the output is wrong from a physical perspective. I am using bvp5c and I'm wondering how to correctly write the residual vector. As my solution vector has 16 components and I have three regions, the solutions passed to bcfun (YL and YR) are 3x16 matrices; therefore the residual vector should have 48 components.

I have 12 boundary conditions; 2 at each outer boundary and 4 at each inner interface. Currently I stack them in the residual vector as follows:

% at -G

res(1:4) = reshape(u1_L + so1, [4, 1]);

res(5:8) = reshape(tu1_L + tso1 , [4, 1]);

% at -L/2

res(9:12) = reshape(gamma1_R - gamma2_L, [4, 1]);

res(13:16) = reshape(tgamma1_R - tgamma2_L, [4, 1]);

res(17:20) = reshape(u1_R - u2_L + so2 , [4, 1]);

res(21:24) = reshape(tu1_R - tu2_L + tso2, [4, 1]);

% at L/2

res(25:28) = reshape(gamma2_R - gamma3_L, [4, 1]);

res(29:32) = reshape(tgamma2_R - tgamma3_L, [4, 1]);

res(33:36) = reshape(u3_L - u2_R + so3, [4, 1]);

res(37:40) = reshape(tu3_L - tu2_R + tso3, [4, 1]);

% at G

res(41:44) = reshape(u3_R + so4, [4, 1]);

res(45:48) = reshape(tu3_R + tso4, [4, 1]);

where gamma and tgamma are the functions I'm looking for and u, tu are their derivatives respectively. So my question is if the order matters at all; I assumed I should go from left to right in the x coordinate and also first write the functions bc and then derivatives bc. To reiterate, is there a general rule when writing down the residual vector?

[Torsten]

We need a mathematical description of the problem (equations, boundary conditions) and the complete code in order to help.

Did you study this page and transfer it to your problem:

https://uk.mathworks.com/help/matlab/math/solve-bvp-with-multiple-boundary-conditions.html

?

If you have 4 first-order differential equations and 3 regions, the matrices YL and YR will be of size 4x3.

[kristjan]

My equations are the following:

with boundary conditions:

1. at outer boundaries:
2. at inner interfaces: ,

where L/R denote outer regions and M denotes the middle one. The differential equation also change from region to region since A = 0 in the middle and the S subscripted quantities are different in each region.

In the equations above I'm solving for and their derivatives. The matrices also have explicit dependence on . All of them are 4x4 matrices so I think that I should have 16x3 dimensions of YL and YR.

I have also studied the webpage and while I understand the example, it's also a quite simple problem. Conceptually my code flows in the same way, just with the solution being a vector of 16 components which are then reshaped into matrices for the purposes of computation in odefun and bcfun. My main concern is just how exactly to stack the residual vector; does it even matter?

function res = bcfun(yL, yR, params)
    % yL represents all solutions on the left of each region and yR
    % conversely on the right side of each region   

    A_x = params.A_x;
    A_cx = params.A_cx;    

    % Extract gamma1 and tilde gamma1 on both boundaries of region 1 (-G &
    % -L/2)
    gamma1_L = reshape(yL(1:4, 1), [2, 2]);         % gamma1 on left boundary of region 1
    tgamma1_L = reshape(yL(5:8, 1), [2, 2]);     % gamma1 tilde on left boundary of region 1
    u1_L = reshape(yL(9:12, 1), [2, 2]);            % gamma1' on left boundary of region 1
    tu1_L = reshape(yL(13:16, 1), [2, 2]);       % til_gamma1' on left of region1

    gamma1_R = reshape(yR(1:4, 1), [2, 2]);         % gamma1 on right boundary of region 1
    tgamma1_R = reshape(yR(5:8, 1), [2, 2]);     % gamma1 tilde on right boundary of region 1
    u1_R = reshape(yR(9:12, 1), [2, 2]);            % gamma1' on right boundary of region 1
    tu1_R = reshape(yR(13:16, 1), [2, 2]);       % tilde gamma1' on right boundary of region 1

    % extract gamma2, tilde gamma2 and their derivatives both boundaries
    % (-L/2 & L/2) 
    gamma2_L = reshape(yL(1:4, 2), [2, 2]);         % gamma2 on the left boundary of region 2
    tgamma2_L = reshape(yL(5:8, 2), [2, 2]);     % tilde gamma2 on the left boundary of region 2
    u2_L = reshape(yL(9:12, 2), [2, 2]);            % gamma2' on left boundary of region 2
    tu2_L = reshape(yL(13:16, 2), [2, 2]);       % tilde gamma2' on left boundary of region 2

    gamma2_R = reshape(yR(1:4, 2), [2, 2]);         % gamma2 on the right boundary of region 2
    tgamma2_R = reshape(yR(5:8, 2), [2, 2]);     % tilde gamma2 on the right boundary of region 2
    u2_R = reshape(yR(9:12, 2), [2, 2]);            % gamma2' on right boundary of region 2
    tu2_R = reshape(yR(13:16, 2), [2, 2]);       % tilde gamma2' on right boundary of region 2
    
    % Extract gamma3 and tilde gamma3 on both boundaries of region 3 (L/2 &
    % G)
    gamma3_L = reshape(yL(1:4, 3), [2, 2]);         % gamma3 on left boundary of region 3
    tgamma3_L = reshape(yL(5:8, 3), [2, 2]);     % gamma3 tilde on left boundary of region 3
    u3_L = reshape(yL(9:12, 3), [2, 2]);            % gamma3' on left boundary of region 3
    tu3_L = reshape(yL(13:16, 3), [2, 2]);       % tilde gamma3' on left boundary of region 3

    gamma3_R = reshape(yR(1:4, 3), [2, 2]);         % gamma3 on right boundary of region 3
    tgamma3_R = reshape(yR(5:8, 3), [2, 2]);     % gamma3 tilde on right boundary of region 3
    u3_R = reshape(yR(9:12, 3), [2, 2]);            % gamma3' on right boundary of region 3
    tu3_R = reshape(yR(13:16, 3), [2, 2]);       % tilde gamma3' on right boundary of region 3  
    
    % Boundary conditions
    
    % define SO terms here so they can be commented out easier
    so1 = - 1i*(A_x*gamma1_L + gamma1_L*A_cx);
    tso1 = 1i*(A_cx*tgamma1_L + tgamma1_L*A_x);
    so2 = - 1i*(A_x*gamma1_R + gamma1_R*A_cx);
    tso2 = + 1i*(A_cx*tgamma1_R + tgamma1_R*A_x);
    so3 = - 1i*(A_x*gamma3_L + gamma3_L*A_cx);
    tso3 = + 1i*(A_cx*tgamma3_L + tgamma3_L*A_x);
    so4 = - 1i*(A_x*gamma3_R + gamma3_R*A_cx);
    tso4 = 1i*(A_cx*tgamma3_R + tgamma3_R*A_x);
    
    % at -G
    res(1:4) = reshape(u1_L + so1, [4, 1]);
    res(5:8) = reshape(tu1_L + tso1 , [4, 1]);
    % at -L/2
    res(9:12) = reshape(gamma1_R - gamma2_L, [4, 1]);
    res(13:16) = reshape(tgamma1_R - tgamma2_L, [4, 1]);    
    res(17:20) = reshape(u1_R - u2_L + so2 , [4, 1]);    
    res(21:24) = reshape(tu1_R - tu2_L + tso2, [4, 1]);
    % at L/2
    res(25:28) = reshape(gamma2_R - gamma3_L, [4, 1]);
    res(29:32) = reshape(tgamma2_R - tgamma3_L, [4, 1]);    
    res(33:36) = reshape(u3_L - u2_R + so3, [4, 1]);
    res(37:40) = reshape(tu3_L - tu2_R + tso3, [4, 1]);
    % at G    
    res(41:44) = reshape(u3_R + so4, [4, 1]);
    res(45:48) = reshape(tu3_R + tso4, [4, 1]);  
    
end