Draw this fractal generated by applying Newton's method to cosh(x) - 1

Question

I came across this picture the other day: (Credit to Josep M Batlle I Ferrer)

Your job is to generate this picture. This graph is generated by repeatedly applying newton's method to the graph of:

$$f(x)=\cosh(x)-1$$

In other words, repeatedly apply the following function to every point till the rectangular distance is less than some limit to a attractor:

$$g(x)=\frac{x \cdot \sinh(x) - \cosh(x) + 1}{\sinh(x)}$$

The process is described in psuedocode below:

let n be the maximum number of iterations
for every point (as a complex number) in the image:
    repeat n times:
        point = (point * sinh(point) - cosh(point) + 1) / sinh(point)
        
        if the rectangular distance between the point and any attractor is less than 0.01:
             break
        end if
    end repeat
    color the point in the image depending on the nearest attractor and the number of iterations required

The rectangular distance is defined as max(abs(dx), abs(dy)). We use rectangular distance instead of straight line distance so we get the wave shapes that make this graph so beautiful.

There is a attractor at 0+n*2i*PI for any integer N. Thus at even increments of 2 PI in the imaginary direction.

Rules

• Your graph must include the area from [-0.75, -0.75 + PI] to [0.75, +0.75 + PI].
• Resolution over this needs to be at least 256 by 256. Either characters or pixels. If you decide to cover a larger area you need to increase the resolution in proportion.
• Use at least 20 iterations.
• Both adjacent iteration counts and adjacent attractors need to be drawn in distinct colors. You must use at least 3 different colors. Coloring should soley be based on the attractor and iteration count.
• If any point causes a division by 0 or another math error, you may choose to render anything at this point. This is considered undefined behavior.
• Displaying a image, saving a image, or drawing ASCII art are all acceptable, as long as the other rules are met. If you draw ASCII art characters need to be ANSI colored.
• Please include a screenshot of the result in your answer

This is code golf, shortest answer wins. ASCII art answers do not compete against graphical answers in the same language.

Background info

• Newton's method
• Newton's fractal
• Hyperbolic Trignometry
• Julia set

Answer 1

Desmos, 401 375 357 353 bytes

-26 bytes thanks to @aiden-chow

-18 bytes thanks to flexabrotnt#1409

H(u,v)=(u.xv.x-v.yu.y,u.xv.y+v.xu.y)
F(u,v)=(u.xv.x+u.yv.y,u.yv.x-u.xv.y)/(v.x²+v.y²)
E(z)=e^{z.x}(cos(z.y),sin(z.y))
C(z,k)=H((1/2,0),E(z)+kE(H((-1,0),z)))
N(z)=F(H(z,C(z,-1))-C(z,1)+(1,0),C(z,-1))
G(z)=max(z.x²,z.y²)
Y(z)=N(N(N(N(N(z)))))
Z(z)=Y(Y(Y(Y(z))))
l=[0...2]
O=hsv(l110,1,1)
P=Z((x,y))
X=round(P.y/(2π))
G((0,2Xπ)-P)<0.01\{\mod(X,3)=l\}

Try It On Desmos!

Before golfing

Answer 2

Rust + Num ASCII Art - 420 413 bytes

-6 bytes thanks to @LF

use{num::complex::*,std::f32::consts::*};(0..255).any(|y|(0..512).any(|x|if let Some(n)=(0..99).scan(Complex::new(x as f32/256.0-1.,y as f32/256.0-0.5+PI),|s,l|{*s=*s-(s.cosh()-1.)/s.sinh();Some((*s,l))}).filter(|s|s.0.re.abs().max(s.0.im-((s.0.im/(2.*PI)).round()*2.*PI).abs())<0.01).next(){print!("\x1b[{:0<2};{:0<2}m0",30+n.1%8,40+((n.0.im/PI/2. + 36.).round()%8.)as usize)}else{print!("0")}>())||print!("
")>());

Answer 3

Python + Pygame, 410 392 bytes (Example Submission)

-16 bytes thanks to wizzwizz4 and others

from pygame import*
from cmath import*
w=256;s=Surface((w,w));f=lambda p:p-(cosh(p)-1)/sinh(p);g=lambda t:round(t.imag/pi/2)*pi*2j;h=lambda p:max(abs(p.real),abs(p.imag))>0.01
for x in range(w):
 for y in range(w):
  t=x/w-0.5+(y/w-0.5-pi)*1j
  for i in range(99):
   try:t=f(t);1/h(g(t)-t)
   except:break
  s.set_at((x,y),[0,[w-1,5,5],[9,w-1,9]][(round(t.imag/pi/2)+i)%3])
image.save(s,"f")

Result will be saved as a TGA file named "f":

Answer 4

JavaScript (ES6), 271 260 253 251 bytes

c=>{z=c.getContext`2d`;for(i=384;i--;)for(j=384;i--;)with(Math){x=i/256-.75;y=j/256-.75+PI;for(k=10;k<30>(x*x<1e-4&abs(y-2*PI*(v=round(y/2/PI)))<.01);k++)w=cosh(x)+cos(y),x-=sinh(x)/w,y-=sin(y)/w;z.fillStyle=`#0${(k+v%10)*3%10}9`;z.fillRect(i,j,1,1)}}

Takes as parameter a 384×384 canvas on which to draw the result. Demo, with whitespace for readability and animation:

f = async c => {
  z = c.getContext`2d`;
  for (i = 384; i--; ) {
    for (j = 384; j--; ) with (Math) {
      x = i / 256 - .75;
      y = j / 256 - .75 + PI;
      for (k = 10; k < 30 > (x * x < 1e-4 & abs(y - 2 * PI * (v = round(y / 2 / PI))) < .01); k++)
        w = cosh(x) + cos(y),
        x -= sinh(x) / w,
        y -= sin(y) / w;
      z.fillStyle = `#0${(k+v%10)*3%10}9`;
      z.fillRect(i, j, 1, 1);
    }
    await new Promise(resolve => requestAnimationFrame(resolve));
  }
};
f(c);

<canvas width=384 height=384 id=c>

Edit: Saved 11 bytes thanks to @Ausername. Saved 7 bytes thanks to @Steffan. Saved 2 bytes thanks to @Dingus.

Answer 5

Pyxplot 0.8.4, ^{190 188} 186 bytes

se nu c
su g(x,y){z=x+y*i
fo n=1to20{z=z-coth(z)+csch(z)
if(Re(z)**2<1e-4)*cos(Im(z))>cos(.01){ret ceil(Im(z)/pi-.5);}
}
}
se sa gr 1e3x1e3
se colm hsb c2:c1:1
p[-1:1][2:4]g(x,y):n w col

Conveniently for code golf, most commands in Pyxplot may be abbreviated zealously. By default (empty .pyxplotrc config file), output is written to pyxplot.eps.

If the text and axis labels produced by default must be removed, replace the last line with

se ax xy in
p[-1:1][2:4]g(x,y):n t''w col

at a cost of 15 bytes. Add se s s for a square aspect ratio.

Some mathematical simplifications are used. Firstly, we write \$g\$ as \$g(z) = z-\coth(z)+\mathop{\rm csch}(z)\$. Secondly, we note that \$\max(\lvert\Delta x\rvert, \lvert\Delta y \rvert) < t\$ implies that both \$\lvert\Delta x\rvert < t\$ and \$\lvert\Delta y\rvert < t\$ and treat these two constraints separately:

• \$\lvert\Delta x\rvert = \lvert x-0\rvert = \lvert x\rvert\$, so the first constraint becomes \$x^2<t^2\$, which is shorter in Pyxplot than \$\lvert x\rvert < t\$ (which would be abs(Re(z))<.01).
• \$\lvert\Delta y\rvert = \lvert y-2k\pi\rvert\$ for some integer \$k\$. Observe that the attractors lie at points where the cosine function is maximal, i.e. \$\cos(2k\pi)=1\$. Thus the constraint on \$y\$ may be written as \$\cos(\lvert y-2k\pi\rvert) > \cos(t)\$, which simplifies to \$\cos(y) > \cos(t)\$ by the symmetry and periodicity of the cosine function.

For the plot, an HSB colourmap is used. At each point the hue is determined by the number of iterations (n in the code), automatically rescaled to span the range from 0 and 1 (i.e. the full range of hues). The saturation is determined by the value of \$2k\$ (the return value of g(x,y)) for the point's attractor, again rescaled to span the full range of saturations. The brightness is fixed at maximum.

An iteration limit of 20 is imposed. Red and orange hues become more predominant at the expense of green and blue as this limit is increased.

In theory, it is possible to omit [-1:1][2:4] in the last line of the code (saving 10 bytes), have Pyxplot apply its default axis ranges (\$-10\le x,y\le 10\$), and increase the grid resolution (1e3x1e3) to compensate. In practice, setting the resolution to 4e3x4e3 exhausted the available memory (16 GB RAM plus 20 GB swap) on my machine.

Ungolfed code

set numerics complex
subroutine g(x, y) {
  z = x + y * i
  for n = 1 to 20 {
    z = z - coth(z) + csch(z)
    if Re(z)**2 < 1e-4 and cos(Im(z)) > cos(1e-2) {
      return ceil(Im(z) / pi - 0.5)
    }
  }
}
set samples grid 1e3 x 1e3
set colourmap hsb c2:c1:1
plot [-1:1] [2:4] g(x, y):n with colourmap

Answer 6

Wolfram Language (Mathematica), 128 bytes

Image@Table[n=0;x+y*I//.z_/;Max@Abs@{Re@z,Mod[i=Im@z,Pi,n+=8t;-1]}>2t:>z-Coth@z+Csch@z;i~Hue~n,{y,-1+Pi,1+Pi,t=.005},{x,-1,1,t}]