Van der Pol Oscillator (https://github.com/hiroloquy/van-der-pol-oscillator.git / https://youtu.be/KzInKpTw0MY)

vdp_6x.plt

reset

#=================== Setting ====================
set nokey
L = 5.0
set xr[-L:L]
set yr[-L:L]
set xl "{/Times:Italic=18 x}({/Times:Italic=18 t})" font 'Times New Roman, 18'
set yl "{/Times:Italic=18 x}'({/Times:Italic=18 t})" font 'Times New Roman, 18'
set tics font 'Times New Roman,14'
set xtics 1
set ytics 1
set size ratio 1
unset grid

# Output GIF file
set term gif animate delay 4 size 1280,720
set output sprintf("vdp_6x_mu=%0.2f.gif", mu)

#=================== Parameter ====================
mu = 1.         # Coefficient of VDP equation
dt = 0.001      # Time step [s]
dh = dt/6.      # Coefficient of Runge-Kutta 4th
lim1 = 70       # Stop time
end = 30        # Time limit [s]
lim2 = end/dt   # Number of calculation
cut = 100       # Decimation

off = 0.15      # Offset
r = 0.08        # Size of points

N    = 6        # Number of points (array size)
array x[N]      # Position of points
array y[N]

# Color of points
array c[N] = ["red", "royalblue", "spring-green", \
              "dark-pink", "grey80", "dark-orange"]

#=================== Function ====================
#Van der Pol equation
f1(x, y)   = y
f2(x, y)   = mu * (1 - x**2) * y - x

#Runge-Kutta 4th order method
rk4x(x, y) = (k1 = f1(x, y),\
              k2 = f1(x + dt*k1/2, y + dt*k1/2),\
              k3 = f1(x + dt*k2/2, y + dt*k2/2),\
              k4 = f1(x + dt*k3, y + dt*k3),\
              dh * (k1 + 2*k2 + 2*k3 + k4))
rk4y(x, y) = (k1 = f2(x, y),\
              k2 = f2(x + dt*k1/2, y + dt*k1/2),\
              k3 = f2(x + dt*k2/2, y + dt*k2/2),\
              k4 = f2(x + dt*k3, y + dt*k3),\
              dh * (k1 + 2*k2 + 2*k3 + k4))

# Time and parameter(mu)
Label(t, mu) = sprintf("{/Times:Italic t} = %.1f\n{/symbol:Italic m} = %.1f", t, mu)

#=================== Plot ====================
# Initiate value
t  = 0.           # Time
x[1] = 0.5        # x1(t)
y[1] = 0.0        # x1'(t)
x[2] = 1.0        # x2(t)
y[2] = -3.0       # x2'(t)
x[3] = -1.7       # x3(t)
y[3] = 1.9        # x3'(t)
x[4] = 1.0        # x4(t)
y[4] = 2.1        # x4'(t)
x[5] = 3.5        # x5(t)
y[5] = 3.7        # x5'(t)
x[6] = -3.        # x6(t)
y[6] = -3.        # x6'(t)

# Axes
set arrow 1 from -L, 0 to L, 0 lc -1 lw 2 back
set arrow 2 from 0, -L to 0, L lc -1 lw 2 back

# Draw initiate state for lim1 steps
do for [i = 1:lim1] {
    # Display time and mu
    set label 1 left Label(t, mu) font 'Times New Roman, 18' at graph 0.02, 0.95

    # Display coordinates of N points
    if(i < lim1*0.8){
        do for[j = 1:N] {
            set label j+1 left sprintf("(%.1f, %.1f)", x[j], y[j])\
                font 'Times New Roman, 14' at x[j]+off, y[j]+off
        }
    }

    # Hide the coordinates
    if(i == lim1*0.8){
        do for[j = 1:N] {
            unset label j+1
        }
    }

    # N points
    do for[j=1:N] {
        set obj j circ at x[j], y[j] size r fc rgb c[j] fs solid border -1 lw 2 front
    }

    plot 1/0
}

# Update for lim2 steps
do for [i = 1:lim2] {
    t = t + dt

    # Calculate using Runge-Kutta 4th
    do for[j = 1:N]{
        x[j] = x[j] + rk4x(x[j], y[j])
        y[j] = y[j] + rk4y(x[j], y[j])
    }

    # Time and mu
    set label 1 Label(t, mu)

    # Objects
    do for[j = 1:N] {
        set obj j at x[j], y[j]
        set obj 6*i+j circ at x[j], y[j] size 0.09*r fc rgb c[j] fs solid
    }

    # Decimate and display
    if(i%cut==0){
        plot 1/0
    }
}

set out

vdp_equation_gif.plt

reset

#=================== Parameter ====================
mu = 1.0		# Coefficient of VDP equation
dt = 0.001		# Time step [s]
dh = dt/6.		# Coefficient of Runge-Kutta 4th
lim1 = 100		# Stop time
end = 40		# Time limit [s]
lim2 = end/dt	# Number of calculation
cut = 100		# Decimation

#=================== Function ====================
f1(a, b) = b
f2(a, b) = mu * (1 - a**2) * b - a

title(t, mu) = sprintf("{/Times:Italic t} = %.1f\n{/symbol:Italic m} = %.1f", t, mu)

#=================== Setting ====================
set term gif animate delay 4 size 426, 426
filename = sprintf("vdp_mu=%0.2f.gif", mu)
set output filename

set nokey
L = 3.0
set xr[-L:L]
set yr[-L:L]
set xl "{/Times:Italic=18 x}({/Times:Italic=18 t})" font 'Times New Roman, 18'
set yl "{/Times:Italic=18 x}'({/Times:Italic=18 t})" font 'Times New Roman, 18'
set tics font 'Times New Roman,14'
set xtics 1
set ytics 1

set size ratio 1
unset grid

#=================== Plot ====================
# Initial value
x = 0.05		# x(t)
y = 0.			# x'(t)
t = 0.

position = sprintf("(%.2f, %d)", x, y)

# Draw initiate state for lim1 steps
do for [i = 1:lim1] {
    # Display time and mu
	set label 1 right title(t, mu) font 'Times New Roman, 14' at graph 0.98, 0.95

    # Display coordinates of N points
	if(i < lim1*0.8){
		set label 2 left position font 'Times New Roman, 14' at x+0.2, y+0.2
	} else {
		if(i==lim1*0.8){
			unset label 2
		}
	}

	set arrow 1 from -L, 0 to L, 0 lc -1 lw 2 back
	set arrow 2 from 0, -L to 0, L lc -1 lw 2 back
	set object 1 circle at x, y size 0.06 fc rgb 'black' fs solid back
	set object 2 circle at x, y size 0.03 fc rgb 'red' fs solid front

	plot 1/0
}

# Update for lim2 steps
do for [i = 1:lim2] {
	t = t + dt

    # Calculate using Runge-Kutta 4th
	k11 = f1(x, y)
	k12 = f2(x, y)
	k21 = f1(x + dt*k11/2.,	y + dt*k12/2.)
	k22 = f2(x + dt*k11/2.,	y + dt*k12/2.)
	k31 = f1(x + dt*k21/2.,	y + dt*k22/2.)
	k32 = f2(x + dt*k21/2.,	y + dt*k22/2.)
	k41 = f1(x + dt*k31,	y + dt*k32)
	k42 = f2(x + dt*k31,	y + dt*k32)

	x = x + dh * (k11 + 2*k21 + 2*k31 + k41)
	y = y + dh * (k12 + 2*k22 + 2*k32 + k42)

	# Time and mu
	set label 1 title(t, mu)

	# Plot
	set object i+1 size 0.002 behind
	set object 1 at x, y
	set object i+2 circle at x, y size 0.03 fc rgb 'red' fs solid front

	if(i%cut==0){
		plot 1/0
	}
}

set out

vdp_vary_mu.plt

reset

#=================== Parameter ====================
# mu = 1.			# Coefficient of VDP equation
dt = 0.004		# Time step [s]
dh = dt/6.		# Coefficient of Runge-Kutta 4th
lim1 = 60		# Stop time
end = 100		# Time limit [s]
lim2 = end/dt	# Number of calculation
cut = 500		# Decimation

#=================== Function ====================
f1(a, b) = b
f2(a, b) = mu * (1 - a**2) * b - a

title(t, mu) = sprintf("{/Times:Italic t} = %.1f\n{/symbol:Italic m} = %.1f", t, mu)

#=================== Setting ====================
set term gif animate delay 4 size 1280,720
filename = "vdp_mu_10patterns.gif"
set output "van der Pol equation 2.gif"

set nokey
L = 8.0
set xr[-L:L]
set yr[-L:L]
set xl "{/Times:Italic=22 x}({/Times:Italic=22 t})" font 'Times New Roman, 22'
set yl "{/Times:Italic=22 x}'({/Times:Italic=22 t})" font 'Times New Roman, 22'
set tics font 'Times New Roman,18'
set xtics 2
set ytics 2

set size ratio 1
unset grid

#=================== Plot ====================
do for [j=0:50:5]{
	# Change the coefficient of VDP equation
	if(j==0){
		mu = 0.1
	} else {
		mu = j/10.
	}

	# Initial value
	x = 0.1		# x(t)
	y = 0.		# x'(t)
	t = 0.

	position = sprintf("(%.1f, %d)", x, y)

	# Draw initiate state for lim1 steps
	do for [i = 1:lim1] {
		# Display time and mu
		set label 1 right title(t, mu) font 'Times New Roman, 22' at graph 0.98, 0.95
		
    	# Display coordinates of N points
		if(i < lim1*0.8){
			set label 2 left position font 'Times New Roman, 14' at x+0.2, y+0.2
		} else {
			if(i==lim1*0.8){
				unset label 2
			}
		}

		set arrow 1 from -L, 0 to L, 0 lc -1 lw 2 back
		set arrow 2 from 0, -L to 0, L lc -1 lw 2 back
		set object 1 circle at 0, 0 size 2 fc rgb 'black' fs solid behind
		set object 2 circle at 0, 0 size 1.95 fc rgb 'white' fs solid behind
		set object 3 circle at x, y size 0.04 fc rgb 'black' fs solid back
		set object 4 circle at x, y size 0.02 fc rgb 'red' fs solid front

		plot 1/0
	}

	# Update for lim2 steps
	do for [i = 1:lim2] {
		t = t + dt

		# Calculate using Runge-Kutta 4th
		k11 = f1(x, y)
		k12 = f2(x, y)
		k21 = f1(x + dt*k11/2.,	y + dt*k12/2.)
		k22 = f2(x + dt*k11/2.,	y + dt*k12/2.)
		k31 = f1(x + dt*k21/2.,	y + dt*k22/2.)
		k32 = f2(x + dt*k21/2.,	y + dt*k22/2.)
		k41 = f1(x + dt*k31,	y + dt*k32)
		k42 = f2(x + dt*k31,	y + dt*k32)

		x = x + dh * (k11 + 2*k21 + 2*k31 + k41)
		y = y + dh * (k12 + 2*k22 + 2*k32 + k42)

		# Time and mu
		set label 1 title(t, mu)

		# Plot
		set object i+3 size 0.003 behind
		set object 3 circle at x, y
		set object i+4 circle at x, y size 0.02 fc rgb 'red' fs solid front

		if(i%cut==0){
			plot 1/0
		}
	}

	# Pause for lim1 steps
	do for [i=1:lim1] {
		plot 1/0
	}
}
set out