In an earlier post, I described the 1-dimensional Ginzburg-Landau equation and showed how it can be linearized and solved with a implicit differencing scheme. The most interesting feature of the solutions was the appearance of seemingly random oscillations. A similar solution method is possible for the 2d version of the equation:

where again and are real valued constants.

An R language code for solving this with parameter values , and an initial state which is a mixture of 2d plane waves is shown below.

library(graphics) #load the graphics library needed for plotting lx <- 80.0 #length of the computational domain in x-direction ly <- 80.0 #length of the computational domain in y-direction lt <- 60.0 #length of the simulation time interval nx <- 50 #number of discrete lattice points in x-direction ny <- 50 #number of discrete lattice points in y-direction nt <- 240 #number of timesteps dx <- lx/nx #length of one discrete lattice cell in x-direction dy <- ly/ny #length of one discrete lattice cell in y-direction dt <- lt/nt #length of timestep a <- 0 b <- 1.5 kappa1 = dt*(1+1i*a)/dx/dx kappa2 = dt*(1+1i*b) C = c(1:(nx*ny)) Cu = c(1:(nx*ny)) A2d = matrix(nrow=ny,ncol=nx) xaxis <- c(0:(nx-1))*dx #the x values corresponding to the discrete lattice points yaxis <- c(0:(ny-1))*dy #the y values corresponding to the discrete lattice points A = matrix(nrow=(nx*ny),ncol=(nx*ny)) IP = matrix(nrow=4*nx,ncol=4*nx) for (i in c(1:ny)) { for (j in c(1:nx)) { A2d[i,j] <- 0.01*exp(1i*5.21*i+1i*10.331*j)+0.01*exp(1i*15.71*i+1i*17.831*j) } } for (k in c(1:nt)) { #main time stepping loop for(i in c(1:ny)) { for(j in c(1:nx)) { C[(i-1)*nx+j] <- A2d[i,j] } } for(i in c(1:(nx*ny))) { for(j in c(1:(nx*ny))) { A[i,j] <- 0 if(i==j && j!=1 && j!=nx && i!=1 && i!=ny) A[i,j] <- 1+2*kappa1+kappa2*abs(C[j])*abs(C[j]) - dt if(i==j && (j==1 || j==nx) && i!=1 && i!=ny) A[i,j] <- 1+2*kappa1+kappa2*abs(C[j])*abs(C[j]) - dt if(i==j && j!=1 && j!=nx && (i==1 || i==ny)) A[i,j] <- 1+2*kappa1+kappa2*abs(C[j])*abs(C[j]) - dt if(i==j && (j==1 || j==nx) && (i==1 || i==ny)) A[i,j] <- 1+2*kappa1+kappa2*abs(C[j])*abs(C[j]) - dt if(j==i+1 && (i%%nx != 0)) A[i,j] <- -kappa1 if(j==i-1 && (i%%nx != 1)) A[i,j] <- -kappa1 if(j==i+nx) A[i,j] <- -kappa1 if(j==i-nx) A[i,j] <- -kappa1 } } Cu <- solve(A,C) for(i in c(1:ny)) { for(j in c(1:nx)) { if(i==1) Cu[(i-1)*nx+j]=Cu[i*nx+j] if(i==ny) Cu[(i-1)*nx+j]=Cu[(i-2)*nx+j] if(j==1) Cu[(i-1)*nx+j]=Cu[(i-1)*nx+j+1] if(j==nx) Cu[(i-1)*nx+j]=Cu[(i-1)*nx+j-1] } } for(i in c(1:ny)) { for(j in c(1:nx)) { A2d[i,j] <- Cu[(i-1)*nx+j] } } for(l in c(1:(nx-1))) { for(m in c(1:(nx-1))) { #make a bitmap with 4 times more pixels, using linear interpolation IP[4*l-3,4*m-3] = A2d[l,m] IP[4*l-2,4*m-3] = A2d[l,m]+0.25*(A2d[l+1,m]-A2d[l,m]) IP[4*l-1,4*m-3] = A2d[l,m]+0.5*(A2d[l+1,m]-A2d[l,m]) IP[4*l,4*m-3] = A2d[l,m]+0.75*(A2d[l+1,m]-A2d[l,m]) } } for(l in c(1:(4*nx))) { for(m in c(1:(nx-1))) { IP[l,4*m-2] = IP[l,4*m-3]+0.25*(IP[l,4*m+1]-IP[l,4*m-3]) IP[l,4*m-1] = IP[l,4*m-3]+0.5*(IP[l,4*m+1]-IP[l,4*m-3]) IP[l,4*m] = IP[l,4*m-3]+0.75*(IP[l,4*m+1]-IP[l,4*m-3]) } } #make plots of C(x,y) on every third timestep jpeg(file = paste("plot_",k,".jpg",sep="")) image(Re(IP),zlim=c(-3,3)) title(paste("Real part of solution A(x,y,t)",k*dt)) dev.off() }

The code produces 2d plots of the real part of the solution on each timestep, and in the video shown below they have been combined into an animation.

In the animation we see the appearance of spiral patterns typical for these values of parameters . Other values of the parameters

produce different kinds of patterns, as is described in this link.