Too many input arguments while solving BVP using bvp4c

Hello, I am trying to run the following chunk of code, which integrates a system of 11 ODEs with an eigenvalue to be fixed. I get a 'too many arguments' error on calling the function bvp_RHS(). I don't get why:

phif = 0.1; % total length
phit = 0.05;
alpha = 0.3;

C_init = 1e-2;
solinit = bvpinit(linspace(-phif/2, phif/2, 10), @(x)mat4init(x, phif), C_init);

options = bvpset('Stats','on','RelTol',1e-5);
sol = bvp4c(@bvp_RHS, @(yL,yR)bvp_BC(yL, yR, phit, alpha), solinit, options);

%-------------------------------------------------------------------------
function yinit = mat4init(x, phif)
yinit = [cos(pi/phif*x)
         -pi/phif * sin(pi/phif*x)
         0.1
         0.1
         0.1
         -0.1
         0.1
         0
         0
         0
         0
    ];
end

function dyds = bvp_RHS(y,C)
% Total 2+9=11 variables.
dyds = [ y(2)                   
         -0.5*y(1)^3 - C*y(1)   % -0.5*kappa^3 - C*kappa
         y(6)                   % =tx
         y(7)                   % =ty
         y(8)                   % =tz
         -y(3) + y(1)*y(9)      % =-rx + kappa*nx
         -y(4) + y(1)*y(10)     % =-ry + kappa*ny
         -y(5) + y(1)*y(11)     % =-rz + kappa*nz
         -y(1)*y(6)             % =-kappa*tx
         -y(1)*y(7)             % =-kappa*ty
         -y(1)*y(8)             % =-kappa*tz
        ];
end
%-------------------------------------------------------------------------

function res = bvp_BC(yL, yR, phit, alpha)
% BCs at L and R boundaries, with constraint for eigenvalue C.
% Symmetrised BCs for r and t, so that L and R are at \pm phit/2.
res = [ yL(1)
        yR(1)
        yL(3)-cos(phit/2)                 % fix position: rx
        yL(4)+sin(phit/2)                   % fix position: ry
        yL(5)                   % fix position: rz
        yR(3)-cos(phit/2)         % fix position: rx
        yR(4)-sin(phit/2)         % fix position: ry
        yR(5)                   % fix position: rz 
        yL(6)-sin(phit/2)*cos(alpha)                   % fix slope: tx
        yL(7)-cos(phit/2)*cos(alpha)        % fix slope: ty
        yL(8)-sin(alpha)        % fix slope: tz
        yR(6)-sin(phit/2)*cos(alpha)  % fix slope: tx
        ];
end

Could someone please advise?

Answer 1

[Star Strider]

The problem with ‘bvp_RHS’ was actually not with it, but with the extra argument to bvpinit. Eliminating that then produced the error that ‘bvp_BC’ needs to return 11 boundary conditions, not 12 as it currently does. So, I commented-out the last one, and the code ran without error.

Check ‘bvp_BC’ and eliminate the extraneous boundary condition, since I might not have eliminated the correct one (I just wanted to see if the code otherwise worked).

phif = 0.1; % total length
phit = 0.05;
alpha = 0.3;

C_init = 1e-2;
solinit = bvpinit(linspace(-phif/2, phif/2, 100), @(x)mat4init(x, phif));

options = bvpset('Stats','on','RelTol',1e-5);
sol = bvp4c(@bvp_RHS, @(yL,yR)bvp_BC(yL, yR, phit, alpha), solinit, options)

Warning: Unable to meet the tolerance without using more than 909 mesh points.
The last mesh of 892 points and the solution are available in the output argument.
The maximum residual is 0.0191713, while requested accuracy is 1e-05.

The solution was obtained on a mesh of 2674 points.
The maximum residual is  1.917e-02. 
There were 82392 calls to the ODE function. 
There were 435 calls to the BC function. 

sol = struct with fields:
    solver: 'bvp4c'
         x: [1×892 double]
         y: [11×892 double]
        yp: [11×892 double]
     stats: [1×1 struct]

figure
plot(sol.x, sol.y)
grid
legend(compose('$y_{%2d}$',1:11), 'Location','bestoutside', 'Interpreter','latex')

%-------------------------------------------------------------------------
function yinit = mat4init(x, phif)
yinit = [cos(pi/phif*x)
         -pi/phif * sin(pi/phif*x)
         0.1
         0.1
         0.1
         -0.1
         0.1
         0
         0
         0
         0
    ];
end

function dyds = bvp_RHS(C,y)
% Total 2+9=11 variables.
dyds = [ y(2)                   
         -0.5*y(1)^3 - C*y(1)   % -0.5*kappa^3 - C*kappa
         y(6)                   % =tx
         y(7)                   % =ty
         y(8)                   % =tz
         -y(3) + y(1)*y(9)      % =-rx + kappa*nx
         -y(4) + y(1)*y(10)     % =-ry + kappa*ny
         -y(5) + y(1)*y(11)     % =-rz + kappa*nz
         -y(1)*y(6)             % =-kappa*tx
         -y(1)*y(7)             % =-kappa*ty
         -y(1)*y(8)             % =-kappa*tz
        ];
end
%-------------------------------------------------------------------------

function res = bvp_BC(yL, yR, phit, alpha)
% BCs at L and R boundaries, with constraint for eigenvalue C.
% Symmetrised BCs for r and t, so that L and R are at \pm phit/2.
res = [ yL(1)
        yR(1)
        yL(3)-cos(phit/2)                 % fix position: rx
        yL(4)+sin(phit/2)                   % fix position: ry
        yL(5)                   % fix position: rz
        yR(3)-cos(phit/2)         % fix position: rx
        yR(4)-sin(phit/2)         % fix position: ry
        yR(5)                   % fix position: rz 
        yL(6)-sin(phit/2)*cos(alpha)                   % fix slope: tx
        yL(7)-cos(phit/2)*cos(alpha)        % fix slope: ty
        yL(8)-sin(alpha)        % fix slope: tz
%         yR(6)-sin(phit/2)*cos(alpha)  % fix slope: tx
        ];
end

The output doesn’t look very interesting, however if it does what you want, that may not be an issue.

.

[Anshuman Pal]

Thank you very much. I hadn't realised that I did not need an extra BC to determine the eigenvalue C. The output is indeed trivial, but that may be because of the choice of initial guesses. I can work on that.

[Star Strider]

As always, my pleasure!

.