`JumpConditions` needs to be fixed

[trefethen]

@asgeirbirkis, would you take a look at examples/ode/JumpConditions and fix it up, or tell me what to do? Thanks.

[asgeirbirkis]

This is related to Fix explicit jumps for chebops 769.

[tobydriscoll]

Since that chebfun issue is closed, does this example now work?

[asgeirbirkis]

No. One reason is that adchebfun/feval did not support 'left' and 'right' extra arguments. I've added that functionality in 993ae91b on fix-adchebfun-fevalArgs. Still, the following is not working:

N = chebop(@(x,u) 1e-2*diff(u,2) + sin(x).*u);
N.lbc = 1; N.rbc = 1;
N.bc = @(x,u) feval(u,0,'right') - feval(u,0,'left') - 1;
u2 = N\0
plot(u2)
Warning: Matrix is not square! 
> In linop.linsolve at 113
  In chebop.solvebvpLinear at 45
  In chebop.solvebvp at 160
  In chebop.mldivide at 13 
Warning: Matrix is not square! 
> In linop.linsolve at 113
  In chebop.solvebvpLinear at 45
  In chebop.solvebvp at 160
  In chebop.mldivide at 13 
u2 =
   chebfun column (1 smooth piece)
       interval       length   endpoint values  
[      -1,       1]       36         1        1 
Epslevel = 1.330053e-15.  Vscale = 1.669442e+00.

Even when I also explicitly specify the jump on diff(u), a warning is still shown, and the jump does not appear:

>> N.bc = @(x,u) [feval(u,0,'right') - feval(u,0,'left') - 1; ...
    feval(diff(u), 0, 'right') - feval(diff(u), 0, 'left') - 1];
u2 = N\0
Warning: Matrix is not square! 
> In linop.linsolve at 113
  In chebop.solvebvpLinear at 45
  In chebop.solvebvp at 160
  In chebop.mldivide at 13 
Warning: Matrix is not square! 
> In linop.linsolve at 113
  In chebop.solvebvpLinear at 45
  In chebop.solvebvp at 160
  In chebop.mldivide at 13 
u2 =
   chebfun column (1 smooth piece)
       interval       length   endpoint values  
[      -1,       1]       36         1        1 
Epslevel = 1.330053e-15.  Vscale = 1.669442e+00.

[asgeirbirkis]

The following still issues a warning (as expected, since we've agreed not to automatically try to fill in extra jumps conditions where jumps are explicitly imposed):

N.bc = @(x,u) jump(u,0) - 1;
u3 = N\0

But the following is happy:

N.bc = @(x,u) [jump(u,0) - 1; jump(diff(u),0)];
u3 = N\0
u3 =
   chebfun column (2 smooth pieces)
       interval       length   endpoint values  
[      -1,       0]       23         1     0.39 
[       0,       1]       26       1.4        1 
Epslevel = 1.472428e-15.  Vscale = 1.508425e+00.  Total length = 49.

[tobydriscoll]

You can't use feval in this way to impose a jump condition. There's no way for the adchebfun to know that you're taking the difference of right and left in order to update the jump locations. You could hack it to assume a jump location whenever you use either one alone, but that seems to introduce new bug possibilities.

The correct course is to recommend only using jump() to do jump conditions.

[nickhale]

I don't think that would be too much of a hack, but I agree that we should probably just recommend using jump().

[asgeirbirkis]

Got the feval(u,0,'left') syntax working properly in 8e41d9, updating the example now.

[asgeirbirkis]

Sorry, missed out on your comments, I'm happy with just using the jump, and not supporting the feval(u,0,left) syntax.

[tobydriscoll]

Either way, it's a bike shed.

[nickhale]

if you've done it now, we might as well include it and document it as experimental.

[asgeirbirkis]

OK, so have managed to get everything to work up until the last demonstration with multiple jumps:

N = chebop(@(x,u) 1e-2*diff(u,2) + sin(x).*u);
N.lbc = 1; N.rbc = 1;
N.bc = @(x,u) [jump(u,-.8:.2:.8) - (-.8:.2:.8); ...
               jump(diff(u),-.8:.2:.8) + (-.8:.2:.8)];
u4 = N\0
Error using  >= 
Matrix dimensions must agree.
Error in chebDiscretization/whichInterval (line 20)
intNum = find( location >= disc.domain, 1, 'last' );
Error in colloc/feval (line 20)
intnum = disc.whichInterval(location, direction);
Error in functionalBlock>@(z)feval(z,location,direction) (line 211)
            E.stack = @(z) feval(z, location, direction);
Error in functionalBlock>@(z)A.stack(z)+B.stack(z) (line 108)
            C.stack = @(z) A.stack(z) + B.stack(z);
Error in functionalBlock>@(z)A.stack(z)*B.stack(z) (line 93)
                C.stack = @(z) A.stack(z) * B.stack(z);
Error in functionalBlock>@(z)A.stack(z)*B.stack(z) (line 93)
                C.stack = @(z) A.stack(z) * B.stack(z);
Error in colloc/instantiate>instantiateOne (line 46)
    A = disc.source.stack( disc );
Error in colloc/instantiate (line 34)
    M = instantiateOne(disc, data);
Error in colloc/instantiate (line 30)
            M{j,k} = instantiate(discJK);
Error in chebDiscretization/matrix (line 39)
[A, S] = instantiate(disc);
Error in chebDiscretization/getConstraints (line 26)
    B = matrix(disc2, dim, dom);
Error in chebDiscretization/matrix (line 49)
    B = getConstraints(disc);
Error in linop/linsolve (line 110)
        [A, P] = matrix(disc);
Error in chebop/solvebvpLinear (line 45)
u = linsolve(L, rhs, linpref);
Error in chebop/solvebvp (line 160)
    [u, info] = N.solvebvpLinear(L, rhs - residual, pref, displayInfo);
Error in chebop/mldivide (line 13)
[varargout{1:nargout}] = solvebvp(varargin{:}); 

The problem here is that location in chebDiscretization/whichInterval is array-valued. adchebfun seems to be doing the correct thing at least:

>> x = adchebfun('x');
>> u = sin(x);
>> v = jump(u,-.8:.2:.8); v.jumpLocations
ans =
   -0.8000
   -0.6000
   -0.4000
   -0.2000
         0
    0.2000
    0.4000
    0.6000
    0.8000

[asgeirbirkis]

OK, just realised that the last issue is the same as @tobydriscoll mentions at the end of Jump conditions with chebop causing an error, still seems a shame we've lost this capability.