Fractals.h

Go to the documentation of this file.

 
#pragma once
 
#include "lin_alg/newton.h"
 
 
class NewtonFractal
{
public:
  double detail;
  double xMin = -1.0, xMax = 1.0;
  double yMin = -1.0, yMax = 1.0;
  int maxIter = 100;
  double tol = 1e-5;
 
  NewtonFractal(double detailFactor = 1, double _min = -1.0, double _max = 1.0, int _maxIter = 100, double _tol = 1e-5)
    : detail(detailFactor)
    , xMin(_min)
    , xMax(_max)
    , yMin(_min)
    , yMax(_max)
    , maxIter(_maxIter)
    , tol(_tol) { }
 
  Matrix<double> operator()() const {
    size_t n = fabs(xMax - xMin) / (1 / (detail > 10 ? detail : 10));
    auto x_i = linspace(xMin, xMax, n).Transpose();
    auto y_i = linspace(yMin, yMax, n).Transpose();
 
    auto M = zeros(n, n);
    std::vector<Matrix<double>> roots;
 
    for(size_t i = 0; i < n; ++i) {
      for(size_t j = 0; j < n; ++j) {
        auto res = newton(fun, jac, { { x_i(i, 0) }, { y_i(j, 0) } }, tol, maxIter);
        if(res.second == maxIter) {
          M(i, j) = 0;
        } else {
          size_t index = 0;
          for(size_t k = 0; k < roots.size(); ++k) {
            if(norm(res.first - roots[k]) <= tol) { index = k; }
          }
          if(index == 0) {
            index = roots.size();
            roots.push_back(res.first);
          }
          M(i, j) = (double)index;
        }
      }
    }
    return M;
  }
 
private:
  static Matrix<double> fun(const Matrix<double>& x) {
    return Matrix<double>({ { pow(x(0, 0), 3) - 3 * x(0, 0) * x(1, 0) * x(1, 0) - 1 },
                            { -pow(x(1, 0), 3) + 3 * x(0, 0) * x(0, 0) * x(1, 0) } });
  }
  static Matrix<double> jac(const Matrix<double>& x) {
    return Matrix<double>({ { 3 * (x(0, 0) * x(0, 0) - x(1, 0) * x(1, 0)), -6 * x(0, 0) * x(1, 0) },
                            { 6 * x(0, 0) * x(1, 0), 3 * x(0, 0) * x(0, 0) - 3 * x(1, 0) * x(1, 0) } });
  }
};
 
class Mandelbrot
{
public:
  size_t maxIters = 125;
  double startX = -2.5, endX = 1.0;
  double startY = -1.0, endY = 1.0;
  size_t detail = 500;
 
 
  Mandelbrot() { }
 
 
  Matrix<double> operator()() {
    auto stepWidthX = (endX - startX) / (double)detail;
    auto stepWidthY = (endY - startY) / (double)detail;
 
    Matrix<double> M = zeros(detail, detail);
    for(size_t x = 0; x < detail; ++x) {
      for(size_t y = 0; y < detail; ++y) {
        M(x, y) = (double)fun((startX + stepWidthX * (double)x), (startY + stepWidthY * (double)y));
        if(M(x, y) == (double)maxIters) { M(x, y) = 0; }
      }
    }
    return M;
  }
 
private:
  size_t fun(const double& real_c, const double& img_c) {
    size_t iters = 0;
    double real  = 0;
    double img   = 0;
 
    while(++iters < maxIters && real * real + img * img < 4) {
      double temp = real;
      real        = real * real + real_c - img * img;
      img         = 2 * img * temp + img_c;
    }
 
    return iters;
  }
};