How to recreate this result (from a book)?

The result I'm interested in is found within "Synchronization: A Universal Concept in Nonlinear Sciences" page $333$ figure $14.3$. The peculiar fragment is also provided at the end of this post.

So basically there is this dissipative coupling applied to one-dimensional array of initial conditions (horizontal axis) which evolves over time (vertical axis). I realize I won't be able to produce identical result since I don't know the exact initial conditions but that's not the point of this post.

The actual problem is I am not sure how to apply the rule of evolution. If I have some initial conditions which undergo a single iteration the result is essentially a function of those initial conditions... What is this function (procedure) in question?

I'd love to be able calculate it in MatLab. There surely has to be somewhat related standard functions in there...

• You have to solve the given system of PDEs, for $u_1(x,t)$ and $u_2(x,t)$. – nonlinearism Jun 3 '13 at 13:54
• @nonlinearism What about the initial conditions though? I mean some additional information obviously has to be introduced in order to make solution unambiguous. – Pranasas Jun 3 '13 at 19:39
• @WillieWong Thanks for your remark. scicomp it is then! I wasn't quite aware of it. – Pranasas Jun 5 '13 at 15:39
• I've migrated the question over. You may want to re-tag your question. For starters, parabolic-pde and matlab would seem appropriate. – Willie Wong Jun 5 '13 at 15:42
• have you ever tried mathematica...check this demonstrations.wolfram.com/CoupledLorenzOscillators you can specify initial conditions and then function iterates over those. – Rorschach Jun 5 '13 at 18:40

Parabolic PDE's such as those in the book can usually be solved using the Method of Lines. First you create some mesh for the $x$ direciton. I will assume that you used some uniform spacing since the plots don't show any characteristics that show the need of non-uniformity. Next you recast your equations with only the time derviative on the left hand side and change the derivatives in $x$ to finite difference approximations. Here is the first equation for a general interior point:

$\frac{\partial u_{1,i}}{\partial t} = \varepsilon (u_{2,i} - u_{1,i}) - \frac{u_{1,i+1} - 2 u_{1,i} + u_{1,i-1}}{\Delta x^2} - \frac{u_{1,i+2} - 4 u_{1,i+1} + 6 u_{1,i} - 4 u_{1,i-1} + u_{1,i-2}}{\Delta x^4} - u_{1,i} \frac{u_{1,i+1} - u_{1,i-1}}{2 \Delta x}$

I am going to leave the second equation, boundary equations, and points near the boundary to you. Now you have a coupled set of ODE's for $i = 1 ... N$. From your initial conditions you can assign $u_{1,i}$ and $u_{2,i}$ at the first time step. Now at each time step you satisfy the above discretized equations at different times based on which time integration algorithm you are using. If it is explicit Euler, you satisfy it at the beginning of each time step. If it is implicit Euler, the end.

In matlab, however, there is an easy way to handle all of these (and many more complicated) methods. What you want is a function that returns a vector of values equal to the right hand side of the above equation given a vector of $u_{j,i}$. If you assume periodic boundary conditions you get:

function [ u_prime ] = derivative( t, u, delta_x )

u_prime = zeros(length(u),1);
u = [u(end-3:end); u; u(1:4)];

if t < 200
epsilon = 0;
else
epsilon = 0.1;
end

for i = 5:2:length(u) - 4;
u_prime(i-4) = epsilon*(u(i+1) - u(i)) - ...
(u(i+2) - 2*u(i) + u(i-2))/delta_x^2 - ...
(u(i+4) - 4*u(i+2) + 6*u(i) - 4*u(i-2) + u(i-4))/delta_x^4 - ...
u(i)*(u(i+2) - u(i-2))/(2*delta_x);
j = i+1;
u_prime(j-4) = epsilon*(u(j-1) - u(j)) - ...
(u(j+2) - 2*u(j) + u(j-2))/delta_x^2 - ...
(u(j+4) - 4*u(j+2) + 6*u(j) - 4*u(j-2) + u(j-4))/delta_x^4 - ...
u(j)*(u(j+2) - u(j-2))/(2*delta_x);
end
end


Now you can feed this to any of matlabs built-in ODE solvers. I found that ode15s performed ratherly well. I also assumed sinusoidal ICS, but it doesn't appear to matter.

N = 1000;   % Number of space discretizations

x = linspace(0, 150, N);
u_0 = zeros(2*N,1);
u_0(1:2:end-1) = sin(2*x/10);   % u_1
u_0(2:2:end) = -sin(4*x/10);    % u_2
delta_x = x(2) - x(1);

[t,u] = ode15s(@(t,u) derivative(t,u,delta_x), [0 400], u_0);


The results give: