Odd numerical artifact from ode15s where my dependent variable changes despite my derivative vector being zero for that time step.

Hi

I am trying to chase down the origin of this numerical artifact that is messing up my 2d diffusion simulation. It shows up at a corner point in my grid after the first time step. At this point i am approximating the system with two forward 2nd order differences.

if(i==1) 
   ut(i)= D/dx^2*(u(i+2)-2*u(i+1)+u(i)) + ...
              D/dy^2*(u(i+2*mesh_dim(1))-2*u(i+mesh_dim(1))+u(i));

I have a large if elseif structure so the corner point is only calculated once. At the end of pde code I put a break point and checked manually. All the neighboring points involved in the calculation of this point are the same value which would zero the time derivative. When my function is called again at the next time step a change was put in u(1). I needed to change the setting to format long to see it. From there it diffuses through the grid and messes up my solution.

Has anyone else seen this type of behavior and what is the cause. I am trying to step through the solver to identify where it is being introduced but it occurs after ode15s has called my function and was passed back a time derivative vector.

Any ideas?

[SL: edited to correct formatting of text]

Answer 1

[Steven Lord]

"All the neighboring points involved in the calculation of this point are the same value which would zero the time derivative."

How did you confirm this? Did you look at the values as they were displayed or did you use == or subtract the values to confirm they are exactly equal down to the last bit? From the fact that you said:

"I needed to change the setting to format long to see it."

that suggests the neighboring points were not the same value but were very, very close to the same value. Therefore the derivative would be very small but non-zero and the new value should be slightly different from the old value.

[Matthew Rich]

I knew they were exactly the same to the last bit since this is the first time step. I do a global assignment of the initial value of the mesh at the start of my code. I also did the subtraction of the different elements after the last element in my derivative vector is calculated. When I step out of my pde call and I am back in ode15s I see that it has a vector that looks as I would expect. I have traced execution and somewhere around line 125 in the ODE solver. What is confusing is the derivative vector makes sense when the execution leaves my code and comes back messed up.

     if fullJacobian                           % generate full matrix dFdy
  ydel = y(:,ones(1,ny)) + diag(del);
  if isempty(vectvar)
    % non-vectorized
    Fdel = zeros(nF,ny);
    for j = 1:ny
      Fdel(:,j) = feval(F,Fargs{1:diffvar-1},ydel(:,j),Fargs{diffvar+1:end});
    end
    nfcalls = ny;

I see the artifact pop up in the variable fdel. I can post my code if needed.

Thanks in advance.