gradient of a gradient giving drastically different answer than the del2 operator

I took a derivative of a derivative using the gradient method. I then applied the del2 operator to take the second derivative directly. For some reason, the second derivative calculated by taking gradient twice does not at all match that of the del2 operator. The gradient of a gradient method seems to be correct, and the del2 operator just shows a continuously decreasing function (which is wrong). Could anyone take a look at my code and figure out what I am doing wrong with the del2 operator? I also have attached ex.txt so that you can directly compute everything that I am seeing. I also attached the figure that is produced. You can see that the graph on the left (gradient of a gradient) is very different than the del2 method (right graph)

filename = 'ex.txt';
delimiterIn = ' ';
headerlinesIn = 5;
% Import the data file for post-processing
matrix = importdata(filename,delimiterIn,headerlinesIn);
A = matrix.data;
A = A(:,1:3);
%Define the number of distinct isotherms
temp_ids = sum(A(:) == 0.2);
%Define the number of density points sampled
density_ids = length(A)/temp_ids;
%Grab density matrix
density = A((1:density_ids),1);
%Grab temperature matrix
pressure = A(:,3);
%Reshape temperature matrix so it has one column for each temperature
%(isotherm), rather than just a single long column
pressure = reshape(pressure,[density_ids,temp_ids]);
%Differentiate 
%dPdD = gradient(pressure(:,1)) / gradient(density);
[~,dPdD] = gradient(pressure, mean(diff(density)));
[~,ddPddD_grad] = gradient(dPdD, mean(diff(density)));
ddPddD_del2 = 4 *del2(pressure, mean(diff(density)));
subplot(1,2,1)
plot(density, ddPddD_grad)
grid on;
xlim([0.4 0.8])
xlabel('\rho (g/cm^3)');
ylabel('\partial^2p/\partial\rho^2')
subplot(1,2,2)
plot(density, ddPddD_del2)
grid on;
xlim([0.4 0.8])
xlabel('\rho (g/cm^3)');
ylabel('\partial^2p/\partial\rho^2')
temperature = 100:5:300;
density_spacing = density(2,1) - density(1,1); 

Answer 1

[Miriam]

With your gradient method you are only taking the derivative along the y dimension, wheras the del2 operator approximates Laplace's differential operator:

[Benjamin]

Oh I see. Basically, I am trying to take the derivative of pressure with respect to density. I guess the Laplace differential operator does not really give me that. Is there a way to take the second derivative of a function without having to call the gradient operator twice (or diff operator twice)? I am not wantint to differentiate along the x dimension. Each column just represents a different isotherm, and I do not want to differentiate that way. I just want to differentiate pressure with respect to density.

[Miriam]

You can calculate the second derivative using diff once:

ddPddD_diff = diff(pressure,2,1)/mean(diff(density))^2;

although this essentially applies the diff operator twice. I'm not aware of another way of taking second derivatives.

[Benjamin]

Since using the diff operator automatically creates an array that is 1 value less than the function array, is there a way to modify this code to account for that? Currently, I get the error that the vectors must be the same length. Obviously this is because using the diff operator twice makes the array 2 shorter than the array it is being plotted against.

[Miriam]

You can opt to plot the smaller array using:

plot(density(2:end-1), ddPddD_diff);

but if you would like your array of second derivatives to stay the same length, I would just stick with your gradient approach.

[Benjamin]

Great thanks! One last question. Do you know what the difference is here between diff and gradient? They more or less give me the same result, but any "glitches" in the pressure curve get more pronounced using diff than using gradient to calculate the 2nd derivative. Do you know why this is?

[Miriam]

I would check out this answer and the documentation on the two functions.