MATLAB: Solving an ODE for a Complex Variable in Polar Coordinates

complex numbersode45polar coordinates

I want to solve an ODE where the variable of interest is a complex number in terms of polar coordinates, . I have been using the example here where they have chosen . I roughly know what the dynamics of the system should look like (it is a well-known reduction of the Kuramoto model) and am almost certain that it should not have features that cause ode45() to have any problems. However, the systems seems to blow up. I believe it is because of how I have set up definitions or returned values in the Second Function File, but I am not sure. Does anyone have any suggestions? Thanks!
Second Function File and Call:
% Set up ICs and t, then solve
zv0 = [0.8; 3];
tspan = [0 100];
[t, zv] = ode45(@imaginaryODE, tspan, zv0);
% Plot
plot(t,zv(:,1))
function fv = imaginaryODE(t, zv)
% Construct z from argument and angle
z = abs(zv(1)).*exp(1i.*angle(zv(2)));
% Evaluate for the function defined in complexf.m
zp = complexf(t, z);
% Return argument and angle
fv = [abs(zp); angle(zp)];
end
First Function File:
function f = complexf(t, z)
% Constants
omega0 = 2.*pi/24;
gamma = 1/100;
k = 1;
% Function
f = (1i.*omega0 - gamma).*z + (k/2).*(z - (z^2).*conj(z));
end

Best Answer

Hi Cameron,
to solve the differential equation, ode needs real inputs so you could either use z = x+iy in the usual way or use
z = r exp(i psi)            % fwd   
r = abs(z)  psi = angle(z)  % reverse
The problem is with the usual 2pi ambiguity for angles, and Matlab uses -pi < psi <= pi. Angle(z) is discontinuous, which is bad for use in a differential equation.
The following code uses the same x and y approach as in the link you provided (although I changed a bunch of names). For this problem you get nice oscillatory behavior. If you want r and psi, you can always get them later, as in plots 2 and 3. Plot 2 illustrates the problem with angle, and unwrap in plot 3 takes care of it. 
z0 = .8+3i;
xy0 = [real(z0); imag(z0)];
tspan = [0 100];
[t,xy] = ode45(@odefun, tspan, xy0);
figure(1)
x = xy(:,1);
y = xy(:,2);
plot(t,x,t,y); grid on
z = x+i*y;
r = abs(z);
psi = angle(z);
figure(2)
plot(t,r,t,psi); grid on
figure(3)
plot(t,10*r,t,unwrap(psi)); grid on
function Dxy = odefun(t,xy)
% usual ode function,  [x;y] in; [dx/dt; dy/dt] out
z = xy(1) + i*xy(2);            
g = fun(t,z);             
Dxy = [real(g); imag(g)];      
end
function g = fun(t,z)
% Define function that takes and returns complex values
omega0 = 2.*pi/24;
gamma = 1/100;
k = 1;
g = (i*omega0 - gamma)*z + (k/2)*(z - (z.^2).*conj(z));
end
Related Question