BVP with no solution

Consider the equation:

y'' + y = 0

with boundary conditions:

y(0) = 1, y(π) = 1

In this case BVP has no solutions. But when we solve it with bvp4c function we obtain some result:

xlow = 0;
xhigh = pi;
solinit = bvpinit(linspace(xlow,xhigh,20),[0,1]);
sol = bvp4c(@(x,y)[y(2); -y(1)],@(ya,yb)[ya(1); yb(1)-1],solinit);
x = linspace(xlow, xhigh, 100);
y = deval(sol, x);
plot(x, y(1,:))

Moreover, it is not unique:

xlow = 0;
xhigh = pi;
hold all
for k=3:30:200
    solinit = bvpinit(linspace(xlow,xhigh,k),[0,1]);
    sol = bvp4c(@(x,y)[y(2); -y(1)],@(ya,yb)[ya(1); yb(1)-1],solinit);
    x = linspace(xlow, xhigh, 100);
    y = deval(sol, x);
    plot(x, y(1,:))
    drawnow
end

Is there any way to check that solution exists?

Answer 1

[Teja Muppirala]

The boundary conditions you wrote in your code are actually y(0) = 0; y(pi) = 1.

But that's besides the point. As you said yourself, there is no solution either way.

The result from the boundary value solver is garbage. It is very similar to what you get when you try to solve this using the finite difference method:

hold all
for N = 2.^(5:10);
    clear b D
    x = linspace(0,pi,N);
    M = zeros(N,N);
    dx = max(x)/(N-1);
    % Express the equation using a matrix
    for n = 2:(N-1)
        D(n,n+[-1 0 1]) = [1 -2 1]/dx/dx + [0 1 0]; % <-- y'' + y
    end

      % Boundary conditions in the matrix
      D(1,1) = 1;
      D(N,N) = 1;

      b(1) = 0; %x(1) = 0, x(pi) = 1
      b(N) = 1;
      b = b(:);

      %solve for the solution
      y = D\b;
      plot(x,y);
  end

Even though the solution you get actually satisfies the linear equations at every point in the mesh and the boundary conditions are satisfied, as the mesh gets finer and finer, it blows up to infinity (just as you saw in your code when you increased k). This should be a red flag telling you that something is not consistent. In this case, the problem is that you are trying to force out a solution from something that does not have one, and sometimes instead of getting an error, you get erroneous results instead.

Answer 2

[Grzegorz Knor]

I agree with you. Just wonder if you can check the existence of a solution before attempting a solution. The example above is trivial, but often it is very difficult ( in analytical way ) to determine whether the equation has a solution.

My first idea was to check how solution behaves with the change of the grid nodes (k). But it doesn't work perfectly because of the equation with an infinite number of solutions.

For example:

y′′ + y = 0

with boundary conditions:

y(0) = 1, y(2π) = 1

has an infinite number of solutions:

xlow = 0;
xhigh = 2*pi;
hold all
for k=3:30:200
    solinit = bvpinit(linspace(xlow,xhigh,k),[0,1]);
    sol = bvp4c(@(x,y)[y(2); -y(1)],@(ya,yb)[ya(1); yb(1)-1],solinit);
    x = linspace(xlow, xhigh, 100);
    y = deval(sol, x);
    plot(x, y(1,:))
    drawnow
end

But when we are looking for unique solution it seems to be a quite good "start point".

Do you have any other ideas?