How to implement an integral condition when solving a BVP in MatLab

Question

I am trying to solve a coupled system of ODE's using MATLAB's bvp4c function. I want to impose the condition that $$\int_{0}^{\pi} y_{1}(t) y_{1}(t) dt = 1,$$ where $y_{1}$ is the first state variable.

Another post suggested adding a additional equation to the system of ODE's. Namely, $$I(t) = \int_{0}^{t} y_{1}(s) y_{1}(s) ds \, .$$

Then the additional boundary condition $I(\pi)-1=0$ would enforce the integral condition I want. However, I am not sure how to code this in MATLAB. Can anyone help me with this?

As per nicoguaro's comment, the integral equation can be written as $$\frac{d}{dt}I(t) = y_{1}(t) y_{1}(t) \, .$$

Updating the code accordingly, I have

This is the whole code I am using with bvp4c (as kirill said, it would be useful to include the whole code):

s = 100;
solinit = bvpinit(linspace(0,pi,1001),@mat4init,s);
sol = bvp4c(@odeproblem,@eightbc,solinit);

x = linspace(0,pi,1001);
y = deval(sol,x);
trapz(y(1,:).*y(1,:))
y_bar = y(1,:)./norm(y(1,:),2);
plot(x,y(1,:))

fprintf('The intial guess for s was %d. The value solved for was %d.\n',s,sol.parameters)

function dydt = odeproblem(t,y,s)
b = pi;
dydt(1) = y(2);
dydt(2) = y(3);
dydt(3) = y(4);
dydt(4) = y(5);
dydt(5) = y(6);
dydt(6) = y(7);
dydt(7) = y(8);
dydt(8) = (-b^4*y(1) + 2*b^2*y(3) - (b^4-2)*y(5) - 2*s^2*b^2*y(7))/s^2;
dydt(9) = y(1)*y(1);
end

function res = eightbc(ya,yb,s)
res = [ya(1); yb(1); ya(3); yb(3); ya(5); yb(5); ya(7); yb(7); ya(9); yb(9)-1];
end

function yinit = mat4init(x)
yinit = [  cos(2*x), -2*sin(2*x), -4*cos(2*x), 8*sin(2*x),... 
    16*cos(2*x),  -32*sin(2*x), -64*cos(2*x), 128*sin(2*x), 1];
end

Here is a plot of the solution:

All other boundary conditions appear to be met (I didn't include the plots showing that). But the integral condition is still not.

Answer 1

The comments detail a perfectly valid way to solve this by defining the integral as another part of the ODE and adding a boundary condition on that. I'll discuss another way to approach this.

bvp4c is a two-point boundary-value problem solver. Instead of a two-point boundary value problem solver, you need the ability to specify a condition that uses multiple points or, even better, uses a continuous extension of the current solution. You can do something like that using DifferentialEquations.jl because it lets you define residuals with continuous extensions of the solution.

function bc3(residual, sol)
    residual[1] = quadgk((t)->sol(t;idxs=1)^2,0,pi)
    residual[2:8] .= 0 #there's no other boundary conditions mentioned?
end

This works because sol(t) is the solution at time t using an order-matching interpolation.

Okay, so it's possible to handle this directly, but what's the difference? This is the part that's more interesting. Two different cases. When using a shooting-type boundary value problem solver, continuous extensions come naturally and one is doing rootfinding of the initial condition, treating the BVP as an IVP with an unknown initial condition. Thus there is no downside to using this kind of continuous extension condition (provided the BVP solver allows it). But shooting methods of course have a downside that they can't be used for problems which are sensitive to the initial condition, or start with a singularity.

That is not true for collocation or MIRK (Implicit RK) type methods. In these cases, one grids space $x_i$ and relaxes the solution at every point in space $y(x_i)$ via a rootfinder. In this case, the "Jacobian" is the dependency of each point on every other point. Since for Runge-Kutta or collocation methods you only use values from just before to compute the new, this typically has a banded structure. However, that's assuming that the boundaries only depend on $y(x_1)$ and $y(x_{end})$. If you impose a boundary condition which depends on all points in space, like this quadrature form for the integral, then the Jacobian of the full system has a full column which destroys the banded structure.

Since the Jacobian is of size M*N*S where M is the number of spatial points, N is the number of ODEs in the system, and S is the number of stages in the ODE, this can be very large even for simple problems. But as noted before, the two-point BVP problem makes the dependencies such that this is an S-banded matrix, making it easy to solve. bvp4c is one of many such solvers which thus assumes the two-point BVP form because it allows solving this system to be greatly simplified.

If it's not a two-point BVP, this Jacobian is a sparse matrix. Either it's stored actually as a sparse matrix or if the problem is small enough one uses a dense matrix for it, but the inversion algorithm cannot use a banded form and thus resorts to a much slower \. Therefore a solver which generalizes like this should only be used when necessary.

Let me conclude by summarizing this:

1. If you can do a transformation to not have a more complicated boundary condition, that's probably your best bet if it's not difficult to do. Integrals can be added as extra parts in the system of ODEs since an ODE solver is actually just a quadrature method.

2. If the problem is more difficult and such a transformation is not possible, then a shooting-type method is a good bet and doesn't have an overhead for having these kinds of conditions (you just need to find a BVP solver which lets you do it).

3. Shooting methods require that the problem isn't sensitive to the initial condition or have an initial singularity. If that's not the case, you need to use collocation or MIRK-type methods. However, most of these methods, like bvp4c, do not allow this form because they can greatly optimize by assuming the problem has a two-point condition. The few solvers which can handle this case should only be used as a last resort.

Hopefully that explanations helps clear your mind that the workaround from the comments isn't just a hack to get it done, but probably the better way to do it.