MRaster examples 21.0.0.0
Image Processing Library
Loading...
Searching...
No Matches
newton_vs.cpp
Go to the documentation of this file.
1// -*- Mode:C++; Coding:us-ascii-unix; fill-column:158 -*-
2/*******************************************************************************************************************************************************.H.S.**/
3/**
4 @file newton_vs.cpp
5 @author Mitch Richling <https://www.mitchr.me>
6 @brief Draw Fractals via Newton-like methods.@EOL
7 @std C++20
8 @copyright
9 @parblock
10 Copyright (c) 1988-2015, 2017, Mitchell Jay Richling <https://www.mitchr.me> All rights reserved.
11
12 Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
13
14 1. Redistributions of source code must retain the above copyright notice, this list of conditions, and the following disclaimer.
15
16 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions, and the following disclaimer in the documentation
17 and/or other materials provided with the distribution.
18
19 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software
20 without specific prior written permission.
21
22 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
23 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
24 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS
25 OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
26 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH
27 DAMAGE.
28 @endparblock
29 @filedetails
30
31 Laguerre's method:
32 @f[x_k@f]
33 @f[G=\frac{f'(x_n)}{f(x_n)}@f]
34 @f[H=G^2-\frac{f''(x_n)}{f(x_n)}@f]
35 @f[a=\frac{n}{G+-\sqrt{(n-1)(nH-G^2)}}@f]
36 @f[x_{n+1}=x_n-a@f]
37
38 @verbatim
39 f0v=f0(x)
40 if(abs(f0v)<eps) ERROR
41 f1v=f1(x)
42 f2v=f2(x)
43 G=f1v/f0v
44 G2=G*G
45 H=G2-f2v/f0v
46 sqr=(n-1)*(n*H-G2)
47 if(sqr<0) ERROR
48 sq=sqrt(sqr)
49 b1=G+sq
50 b2=G-sq
51 if(abs(b1)>abs(b2))
52 b=b1
53 else
54 b=b2
55 fi
56 if(abs(b)<eps) ERROR
57 a=n/b
58 x=x-a
59 @endverbatim
60
61 Newton's
62 @f[x=x-\frac{f(x_n)}{f'(x_n)}@f]
63
64 @verbatim
65 f0v=f0(x)
66 f1v=f1(x)
67 if(abs(f1v)<eps) ERROR
68 x=x-f0v/f1v
69 @endverbatim
70
71 Maxima:
72 @code{.maxima}
73 f:sin(z);
74 z-f/diff(f, z);
75 @endcode
76
77 Halley's
78 @f[ x_n-\frac{f(x_n)}{f'(x_n)}\left[1-\frac{f(x_n)}{f'(x_n)}\cdot\frac{f''(x_n)}{2f'(x_n)}\right]^{-1} @f]
79
80 @verbatim
81 f0v=f0(x)
82 f1v=f1(x)
83 if(abs(f1v)<eps) ERROR
84 f2v=f2(x)
85 G=f0v/f1v
86 b=1-G*f2v/f1v
87 if(abs(b)<eps) ERROR
88 a=G/b
89 xnext=x-a
90 @endverbatim
91
92*/
93/*******************************************************************************************************************************************************.H.E.**/
94/** @cond exj */
95
96//--------------------------------------------------------------------------------------------------------------------------------------------------------------
97#include "ramCanvas.hpp"
98
99//--------------------------------------------------------------------------------------------------------------------------------------------------------------
100enum class whyStopNV { DIVZERO, //!< Divide by zero (zeroTol)
101 TOOBIG, //!< Iterate got too big (> escapeMod)
102 CONVERGEU, //!< Converged in the upper half plane
103 CONVERGEL, //!< Converged in the lower half plane
104 TOOLONG //!< Too many iterations (> MaxCount)
105 };
106
107enum class solMethNV { NEWTON, //!< Use Newton's method
108 HALLEY, //!< Use Halley's method
109 LAGUERRE //!< Use laguerre's method
110 };
111
112//--------------------------------------------------------------------------------------------------------------------------------------------------------------
113int main(void) {
114 std::chrono::time_point<std::chrono::system_clock> startTime = std::chrono::system_clock::now();
115 const double escapeMod = -32.0;
116 const int MaxCount = 64;
117 const double Tol = 0.0001;
118 const int numToKeep = 1;
119 whyStopNV why;
120 mjr::ramCanvas3c8b theRamCanvas(3840, 2160, -1.5, 1.5, -1.5, 1.5);
121 std::vector<solMethNV> methodsToDo({solMethNV::NEWTON, solMethNV::HALLEY, solMethNV::LAGUERRE});
122
123 for(auto method : methodsToDo) {
124 mjr::ramCanvas3c8b::colorType aColor(255, 255, 255);
125
126 for(int y=0;y<theRamCanvas.getNumPixY();y++) {
127 if(y%512==0)
128 std::cout << "Case: " << (int)method << " Line: " << y << std::endl;
129 for(int x=0;x<theRamCanvas.getNumPixX();x++) {
130 std::complex<double> z = theRamCanvas.int2real(x, y);
131
132 std::vector<std::complex<double>> lastZs(numToKeep);
133 mjr::ramCanvas3c8b::csIntType count = 0;
134 double maxMod = 0.0;
135 while(1) {
136 if (count >= MaxCount) {
137 why = whyStopNV::TOOLONG;
138 break;
139 }
140
141 std::complex<double> f0v=sin(cos(z));
142 std::complex<double> f1v=-sin(z)*cos(cos(z));
143 std::complex<double> f2v=(-sin(z)*sin(z)*sin(cos(z)))-cos(z)*cos(cos(z));
144
145 switch(method) {
146 case solMethNV::NEWTON : {
147 if(abs(f1v) < Tol) {
148 why = whyStopNV::DIVZERO;
149 break;
150 }
151 z=z-f0v/f1v;
152 break;
153 }
154 case solMethNV::HALLEY : {
155 if(abs(f1v) < Tol) {
156 why = whyStopNV::DIVZERO;
157 break;
158 }
159 std::complex<double> G=f0v/f1v;
160 std::complex<double> b=1.0-G*f2v/f1v;
161 if(abs(b) < Tol) {
162 why = whyStopNV::DIVZERO;
163 break;
164 }
165 std::complex<double> a=G/b;
166 z = z-a;
167 break;
168 }
169 case solMethNV::LAGUERRE : {
170 std::complex<double> pdeg = 4.0;
171 if(abs(f0v) < Tol) {
172 break;
173 }
174 std::complex<double> G=f1v/f0v;
175 std::complex<double> G2=G*G;
176 std::complex<double> H=G2-f2v/f0v;
177 std::complex<double> sqr=(pdeg-1.0)*(pdeg*H-G2);
178 std::complex<double> sq=sqrt(sqr);
179 std::complex<double> b1=G+sq;
180 std::complex<double> b2=G-sq;
181 std::complex<double> b=b2;
182 if(abs(b1)>abs(b2))
183 b=b1;
184 if(abs(b) < Tol) {
185 why = whyStopNV::DIVZERO;
186 break;
187 }
188 std::complex<double> a=pdeg/b;
189 z = z-a;
190 break;
191 }
192 }
193
194 double modz = abs(z);
195 lastZs[count%numToKeep] = z;
196 if(modz>maxMod) {
197 maxMod = modz;
198 }
199 if(count >= numToKeep) { // Criteria for convergeance: Last numToKeep z values must be within Tol of each other.
200 bool allEqual = true;
201 for(int i=1; i<numToKeep; i++)
202 if(abs(lastZs[0]-lastZs[i])>Tol) {
203 allEqual = false;
204 break;
205 }
206 if(allEqual) {
207 if(z.real() > 0) {
208 why = whyStopNV::CONVERGEU;
209 } else {
210 why = whyStopNV::CONVERGEL;
211 }
212 break;
213 }
214 }
215 if((escapeMod>0) && (modz>escapeMod)) {
216 why = whyStopNV::TOOBIG;
217 break;
218 }
219 count++;
220 }
221
222 mjr::ramCanvas3c8b::csIntType ccol = (2*4*count);
223 mjr::ramCanvas3c8b::csIntType mcol = ( mjr::ramCanvas3c8b::csIntType )(2*8*maxMod);
224 switch(why) {
225 case whyStopNV::TOOLONG : theRamCanvas.drawPoint(x, y, aColor.cmpRGBcornerDGradiant((mcol)%(2*256), "MWM")); break;
226 case whyStopNV::CONVERGEU : theRamCanvas.drawPoint(x, y, aColor.cmpRGBcornerDGradiant((mcol+ccol)%(2*256), "BWB")); break;
227 case whyStopNV::CONVERGEL : theRamCanvas.drawPoint(x, y, aColor.cmpRGBcornerDGradiant((mcol+ccol)%(2*256), "RWR")); break;
228 case whyStopNV::TOOBIG : theRamCanvas.drawPoint(x, y, mjr::ramCanvas3c8b::colorType(0, 0, 0)); break;
229 case whyStopNV::DIVZERO : theRamCanvas.drawPoint(x, y, aColor.cmpRGBcornerDGradiant((mcol+ccol)%(2*256), "CWC")); break;
230 }
231 }
232 }
233
234 switch(method) {
235 case solMethNV::NEWTON : theRamCanvas.writeTIFFfile("newton_vs_newton.tiff" ); break;
236 case solMethNV::HALLEY : theRamCanvas.writeTIFFfile("newton_vs_halley.tiff" ); break;
237 case solMethNV::LAGUERRE : theRamCanvas.writeTIFFfile("newton_vs_laguerre.tiff"); break;
238 }
239 }
240 std::chrono::duration<double> runTime = std::chrono::system_clock::now() - startTime;
241 std::cout << "Total Runtime " << runTime.count() << " sec" << std::endl;
242}
243/** @endcond */
int main(int argc, char *argv[])