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 '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 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);
y_bar = y(1,:)./norm(y(1,:),2);

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);

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];

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];

Here is a plot of the solution:

Solution from bvp4c

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

  • $\begingroup$ You need to rewrite your new equation as a differential equation as well, an provide it in the same fashion. $\endgroup$
    – nicoguaro
    Apr 19 '17 at 16:43
  • $\begingroup$ Thanks! I think I have correctly done what you said (see edited question), but I am not getting solutions that meet the integral condition. Any suggestions? $\endgroup$
    – BoiseID
    Apr 19 '17 at 21:47
  • $\begingroup$ Can you please include a plot of the solution to make sure nothing weird is going on? You are only checking one of the BCs, but are the other boundary conditions satisfied correctly? Can you add more details, please, including the full invocation of bvp4c? Why do you have 10 BCs (in a function called eightbc) for a TPBVP of size 9? $\endgroup$
    – Kirill
    Apr 19 '17 at 22:50
  • $\begingroup$ $y(t)=[0,0,0,0,0,0,0,0,1]$ is a perfectly valid solution of the BVP that you wrote down, and I think that's the one that Matlab finds, so I think it's not wrong here. $\endgroup$
    – Kirill
    Apr 19 '17 at 23:54
  • $\begingroup$ I've included the full code and a plot of the solution. The BVP requires an additional boundary because of the unknown parameter s. Calling it eightbc was a poor choice of name. It could be that MATLAB is converging to the solution you suggest, but I thought the last equation and BC would prevent that from happening. $\endgroup$
    – BoiseID
    Apr 20 '17 at 1:26

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?

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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.