4

I would like to compute the derivative of a complex-valued function (Holomorphic function) numerically in MATLAB.

I have computed the function in a grid on the complex plane, and I've tried to compute the derivative using the Cauchy–Riemann relations.

Given: u = real(f), v = imag(f), x = real(point), y = imag(point)

The derivative should be given by: f' = du/dx + i dv/dx = dv/dy - i du/dy

where 'd' is the derivative operator.

I've tried the following code:

stepx = 0.01;
stepy = 0.01;
Nx = 2/stepx +1;
Ny = 2/stepy +1;
[re,im] = meshgrid([-1:stepx:1], [-1:stepy:1]);
cplx = re + 1i*im;
z = cplx.^3;

The derivative should be given by: 

f1 = diff(real(z),1,2)/stepx +1i* diff(imag(z),1,2)/stepx;

or

f2 = diff(imag(z),1,1)/stepy - 1i* diff(real(z),1,1)/stepy;

But the two derivatives, which are suppose to be equal, do not match.

What am I doing wrong?

Let's count the number of elements which differs for less than stepx (assuming stepx = stepy):

lm = min(size(f1));
A = f1(1:lm,1:lm);
B = f2(1:lm,1:lm);

sum(sum(abs(A - B) <= stepx))

and using the fix proposed by @A. Donda

f1i = interp1(1 : Ny, f1, 1.5 : Ny);
f2i = interp1(1 : Nx, f2 .', 1.5 : Nx) .';

sum(sum(abs(f1i - f2i) <= stepx))

In the second case they all differ for less than stepx as it should be, while in the first case it's not true.

Simone Gaiarin
  • 340
  • 2
  • 15

1 Answers1

1

The problem is that using the two different expressions, discretized for use in Matlab, you are computing the approximate derivatives at different points in the complex plane. Let's say you are at imaginary value y and you compute the differences along the real axis x, then the ith difference estimates the derivative at (x(i) + x(i + 1))/2, i.e. at all the midpoints between two subsequent x-values. The other way around you estimate the derivative at a given x but at all the midpoints between two subsequent y-values.

This also leads to different sizes of the resulting matrices. Using the first formula you get a matrix of size 201x200, the other of size 200x201. That's because in the first variant there are 200 midpoints along x, but 201 y-values, and vice versa.

So the answer is, you are not doing anything wrong, you are just interpreting the result wrongly.

You can solve the problem by interpolating explicitly along the other dimension (the one not used for the derivative):

f1i = interp1(1 : Ny, f1, 1.5 : Ny);
f2i = interp1(1 : Nx, f2 .', 1.5 : Nx) .';

Where f1 is computed according to your first formula and f2 according to the second. Now both derivatives are evaluated at points that are midpoints along both dimensions, which is why both matrices are of size 200x200.

If you compare them now, you will see that they are identical up to numerical error (after all, diff computes only approximate derivatives and interp1 makes interpolation errors). For your stepsize, this error is maximally 1e-4, and it can be further reduced by using a smaller stepsize.

A. Donda
  • 8,381
  • 2
  • 20
  • 49