title: Why is the total energy of my PDE model not conserved over time?
body: <p>I'm evolving a PDE u(x, y, t) in time using the partial differential equation toolbox in MATLAB <a href="https://i.stack.imgur.com/pLVm9.png" rel="nofollow noreferrer">total energy on y axis vs time on x axis</a>and checking to make sure the total energy in the model is conserved. The energy calculation and equations of motion are found in this paper
  <a href="https://journals.aps.org/prd/pdf/10.1103/PhysRevD.49.2978" rel="nofollow noreferrer">https://journals.aps.org/prd/pdf/10.1103/PhysRevD.49.2978</a>
  What is going wrong in the model that causes total energy to jump around so much?</p>
  <pre><code> b = 1; %constant for shape of potential well
  lambda = 1; %constant for shape of potential well
  width = 4; %adjusts width of initial condition
  sidelength = 30; %sidelength of square geometry
  u0 = @(location) -tanh((location.x.^2 + location.y.^2).^(1./2) - width); %function giving initial conditions
  %u0 = @(location) 2.*exp(-(location.x.^2 + location.y.^2)/width.^2) - 1;
  ut0 = 0; %initial condition of first derivative
  syms u(x, y, t); %field configuration
  syms V(p); %potential equation
  V(p) = (lambda./4) * (p.^2 - b.^2/lambda).^2;
  
  
  dim = 2; %dimensions of space 
  
  %defining time mesh values
  tstart = 0; 
  tstep = .5;
  tend = 50;
  tlist = tstart:tstep:tend; %time mesh
  tPoints = length(tlist);
  
  %https://www.mathworks.com/help/pde/ug/pde.pdemodel.specifycoefficients.html
  %creating 2-dim ODE. To make 3-dim, add z term to u
  ode = diff(u, t, 2) - laplacian(u, [x, y]) == -lambda * u * (u.^2 - (b.^2)/(lambda));
  coeffs = pdeCoefficients(ode, u); %formats pde according to pde toolbox format
  
  
  %https://www.mathworks.com/help/pde/geometry-and-mesh.html
  model = createpde; %creating pde model
  R1 = [3,4,-sidelength/2,sidelength/2,sidelength/2,-sidelength/2,sidelength/2,sidelength/2,-sidelength/2,-sidelength/2]'; %defining geometry by edges 
  gm = [R1];
  g = decsg(gm);
  geometryFromEdges(model,g);
  
  
  
  pdegplot(model, 'EdgeLabels', 'On'); %showing geometry
  
  %https://www.mathworks.com/help/pde/ug/pde.pdemodel.html
  specifyCoefficients(model, &quot;m&quot;, coeffs.m, &quot;d&quot;, coeffs.d, &quot;c&quot;, coeffs.c, &quot;a&quot;, coeffs.a, &quot;f&quot;, coeffs.f);
  %gives correct format of pde to the pde model
  
  %https://www.mathworks.com/help/pde/ug/pde.pdemodel.setinitialconditions.html
  setInitialConditions(model,  u0, ut0); 
  %gives initial conditions to pde model
  
  mesh = generateMesh(model); %automaticcaly meshes geometry 
  pdemesh(model) %shows mesh
  
  
  results = solvepde(model, tlist); %solves pde
  
  sol = results.NodalSolution; %gives pde solutions at each point in time mesh
  
  
  
  
  %-------------------------------Energy calculations---------------------------------------
  
  %mesh(X, Y, uinterp) check to make sure uinterp is correct format
  
  %https://www.mathworks.com/help/pde/ug/pde.stationaryresults.evaluategradient.html#d123e46221
  
  Energy = zeros([1 (tPoints - 2)]); %will contain energies for t= tstart+1 to t = tend - 1
  surfaceEnergies = zeros([1 (tPoints - 2)]);
  kineticEnergies = zeros([1 (tPoints -2)]);
  volumeEnergies = zeros([1 (tPoints -2)]);
  
  numberElements = length(mesh.Elements);
  querypoints = mesh.Nodes;
  mesh = generateMesh(model);
  allElementNodes = (mesh.Elements);
  for t = 2: tPoints - 1
      allValues = sol(:, t); %1 = time
      allValuesn1 = sol(:, t-1);
      allValues1 = sol(:, t+1);
      [gradx, grady] = evaluateGradient(results, querypoints, t);
      
      volumeEnergy = 0;
      kineticEnergy = 0;
      surfaceEnergy = 0;
      
      for elem = 1:numberElements
          elementArea = area(mesh, elem);
          oneElementNodeIndices = allElementNodes(:, elem); 
          oneElementNodeValues = allValues(oneElementNodeIndices);
          oneElementNodeValues1 = allValues1(oneElementNodeIndices);
          oneElementNodeValuesn1 = allValuesn1(oneElementNodeIndices);
          elemValue = sum(oneElementNodeValues) ./ 6;
          elemValue1 = sum(oneElementNodeValues1) ./ 6;
          elemValuen1 = sum(oneElementNodeValuesn1) ./6;
          
          volumeEnergy = volumeEnergy + V(elemValue) .* elementArea;
          
          tdiff = (elemValue1 - elemValuen1) ./ (2*tstep);
          kineticEnergy = kineticEnergy + .5 .* tdiff.^2 .* elementArea;
          
          elemGradx = gradx(oneElementNodeIndices);
          elemGrady = grady(oneElementNodeIndices);
          avgGradx = sum(elemGradx) ./ 6;
          avgGrady = sum(elemGrady) ./6;
          surfaceEnergy = surfaceEnergy + .5 .* (avgGradx.^2 + avgGrady.^2) .* elementArea;
      end
      surfaceEnergies(t-1) = surfaceEnergy;
      kineticEnergies(t-1) = kineticEnergy;
      volumeEnergies(t-1) = volumeEnergy;
      Energy(t-1) = surfaceEnergy + volumeEnergy + kineticEnergy;
      double(Energy(t-1))
  end
  %https://www.mathworks.com/help/pde/ug/pde.femesh-properties.html
  %https://www.mathworks.com/help/pde/ug/pde.pdemodel.generatemesh.html
  
  plot([2:(tPoints - 1)], Energy)
  plot([2:(tPoints - 1)], surfaceEnergies)
  plot([2:(tPoints - 1)], volumeEnergies)
  plot([2:(tPoints - 1)], kineticEnergies)
  save energytest 
  </code></pre>