## Comparison of three quadrature rules

Sometimes, in science and engineering, we need to calculate some definite integrals of the type:

It is equivalent to evaluating the area enclosed between the function graph and x-axis, within the interval [a, b]. There are some integrals (like elliptical, Fresnel, etc.) that cannot be expressed in analytical form. In such cases, we have to use numerical methods, also called quadrature rules. They approximate the integral by a finite sum like follows:

where xi are the abscissas of the nodes, where the function is calculated, and wi are the respective weights. Most quadrature rules use the same formula, but with different values for abscissas and weights.

Probably, most engineers are already familiar with Newton’s or Gauss formulas for numerical integration. However, there are other contemporary methods that can be more efficient for this purpose.

In this post, we will compare three quadrature rules for numerical integration. They will be applied to the same function, and with the same number of nodes. The calculations will be performed with Calcpad, as follows:

Test function (half-circle):

Theoretical value of the integral:

Number of nodes: n = 5

## Boole’s rule

This is one of the Newton-Cotes‘ family of quadrature formulas, with order of n = 4. Nodes are located at equal spacing h, and weights are calculated by approximating the function with Lagrange polynomials.

Spacing – h = 0.5

Integral:

Error:

Scheme:

## Gauss-Legendre’s formula

This is the most simple Gaussian quadrature formula. The function is approximated by Legendre polynomials in the interval of [-1; 1]. Abscissas are determined as the roots of the respective polynomials, so the nodes are not equally spaced.

Abscissas and weights are calculated by the following equations:

Integral:

IG = 2·w1·y1 + 2·w2·y2 + w3·y3 = 2·0.23693·0.42289 + 2·0.47863·0.84265 + 0.56889·1
IG = 1.5759

Error:

Scheme:

This is one of the double exponential integration formulas, proposed by Takahasi and Mori (1974).

Boundary of the interval – ta = 1

Step:

Parameter – tk = –ta + (k – 1) ·h

Formula for calculation of abscissas:

Formula for calculation of weights:

Integral:

Value:

IDE = 2·w1·y1 + 2·w2·y2 + w3·y3 = 2·0.11501·0.30806 + 2·0.48299·0.73848 + 0.7854·1
IDE = 1.5696

Error:

Scheme:

## Summary of the results

The relative errors, from the above three quadrature rules are presented in the following chart:

We can see that the Gaussian quadrature rule is much more precise than the Newton-Cotes’ one of the same order, which is expected. However, Tanh-Sinh quadrature shows the best results from all three. And this is not by chance. According to Wikipedia, this is probably the most efficient numerical integration rule ever known.

That is why, it’s adaptive version is implemented in Calcpad as the main numerical integration method (\$Integral command). Alternatively, you can use the adaptive Gauss-Lobatto quadrature (\$Area command). It is much slower, but unlike the Tanh-Sinh one, it works for functions that are not smooth and continuous.

For the above example, both of them provides extremely high accuracy, that is comparable to the floating point precision of the double (float64) data type:

## How to add title, labels and legend to a Calcpad plot?

Calcpad does not provide options to add titles, axis labels and legends to function plots like in Excel. However, you can decorate your plots with these attributes by your own, by using a little Html.

Bellow, you can find a sample code that you can use as a boilerplate. Since chart colors are predefined, they will be the same for all your future plots:

``````'<div style="width:400pt;">
'<h3><center>Title</center></h3>
'<p>Label Y
\$Plot{f_1(x) & f_2(x) & f_3(x) & f_4(x) & f_5(x) & f_6(x) & x_7|y_7 & x_8|y_8 & x_9|y_9 & x_10|y_10 @ x = a : b}
'<span style="float:right">Label X</span></p>
'<p><b>Legend:</b></p>
'<p><b style="color:Red">▬▬</b> Chart 1 &emsp;&emsp;
'<b style="color:Green">▬▬</b> Chart 2 &emsp;&emsp;
'<b style="color:Blue">▬▬</b> Chart 3</p>
'<p><b style="color:Goldenrod">▬▬</b> Chart 4 &emsp;&emsp;
'<b style="color:Magenta">▬▬</b> Chart 5 &emsp;&emsp;
'<b style="color:DarkCyan">▬▬</b> Chart 6</p>
'<p><b style="color:Purple">●</b> Point 7 &emsp;
'<b style="color:DarkOrange">●</b> Point 8 &emsp;
'<b style="color:Maroon">●</b> Point 9 &emsp;
'<b style="color:YellowGreen">●</b> Point 10</p>
'</div>
``````

And here is how it looks inside Calcpad:

If you have less charts, you can delete the extra lines. Also, there is an option for points, if both x and y are fixed. If you have any questions, please do not hesitate to ask them right away.

## Units in empirical formulas

There are many formulas in structural design codes, that are empirical and if you enter them in Calcpad “as are”, you will not obtain the results in the expected units. Moreover, there are even formulas, that are not dimensionally correct and won’t calculate. Lets take for example the shear resistance without reinforcement, according to Eurocode 2. It is calculated by the following formula:

VRd,c = (CRd,c k (100 ρl   fck)1/3 + k1σcp) bwd (6.2a)

Here, fck and σcp are in MPa, so, we will end up with MPa1/3 + MPa, which does not have physical meaning and cannot be calculated. The k factor is equal to:

k = 1 + sqrt(200/d) ≤ 2.0

If d is in mm, we will get: 1 + mm-0.5, which also does not make any sense.

On the other hand, the minimum resistance is:

VRd,c = (vmin + k1σcp) bwd (6.2b), where

vmin = 0.035k3/2 fck1/2

How can we solve this problem?

Obviously, we have no choice, except to manually compensate for the missing dimensions. For example the value of “200” in the formula for k is actually in mm. So, we will have k = 1 + sqrt(200mm/d) and get a dimensionless result as expected.

If we apply the same approach for the other formulas, we will obtain the following little program for calculation of shear resistance without reinforcement to Eurocode 2:

``````'Section width -'b_w = 250mm
'Effective height -'d = 450mm
'Main reinforcement ratio -'ρ_l = 0.02
'Characteristic concrete strength -'f_ck = 25MPa
'Partial safety factor for concrete -'γ_c = 1.5
'Design concrete strength -'f_cd = f_ck/γ_c
'Normal stress -'σ_cp = 0MPa
'Factors:'k_1 = 0.15', 'C_Rd,c = 0.18/γ_c
k = min(1 + sqr(200mm/d); 2)
v_min = 0.035*k^(3/2)*sqr(f_ck)*MPa^0.5|MPa
'Shear resistance
'<p class="ref">(6.2a)</p>
V_Rd,c = (C_Rd,c*k*(100*ρ_l*f_ck)^(1/3)*MPa^(2/3) + k_1*σ_cp)*b_w*d
'Minimum shear resistance
'<p class="ref">(6.2b)</p>
V_Rd,c = (v_min + k_1*σ_cp)*b_w*d``````

When you run the above program in Calcpad, the results will look as follows:

You should not worry about that formulas do not look exactly as in the design codes. It is far more important to give the correct results as it can be easily verified.

## What’s new in Calcpad 5.8 and later?

In the last several months, Calcpad was subjected to continuous improvement, driven by the efforts of its developers and the help of the community. In this post, we will make a short review of the most important changes.

## 1. Modules

If you have a simple and short worksheet, it is definitely easier to keep it all in a single file. But imagine that you have to build an entire library for design of timber structures. And each worksheet should start with identical content for selection of materials and cross sections. Then, it would be better to put this part in a separate module and reference it from all other worksheets. Now, you can do this in Calcpad, using the following code:

``#include filename.cpd``

If the file is in the same folder, you can enter only the filename, otherwise you must specify the full path. This is better than copying one and the same text in all worksheets. You can also select which part of the module is local and global by inserting the respective keywords: “`#local`“, “`#global`“. Local parts will not be included into the external module.

## 2. String variables and macros

This is a good way to further organize your code and avoid repetitions at a module level. You can assign a piece of code to a string variable. Then, you can insert it wherever you like, by simply writing the name instead of the whole text.

The difference between variables and macros is that macros have parameters, and variables don’t. On the other hand, both can be inline of multiline. All names must and with the “\$” symbol and can include letters, numbers and underscore “_”. You can use the following syntax:

Inline string variable:

``#def variable_name\$ = content``

Multiline string variable:

``````#def variable_name\$
'content line 1
'content line 2
'...
#end def``````

Inline macro:

``#def macro_name\$(param1\$; param2\$; ...) = content``

Multiline macro:

``````#def macro_name\$(param1\$; param2\$; ...)
'content line 1
'content line 2
'...
#end def``````

This is a really powerful feature that can save a lot of work. It allows you to expand the Calcpad functionality with features that never existed before. For example, you can create Html “shortcodes” for faster and better formatting:

``````#def b\$(x\$) = <strong>x\$</strong>
#def i\$(x\$) = <em>x\$</em>
#def sup\$(x\$) = <sup>x\$</sup>
#def sub\$(x\$) = <sub>x\$</sub>
#def err\$(x\$) = <span class="err">x\$</span>
'Usage:
'<p>b\$(Check): i\$(f)sub\$(yk) = err\$(2) kN/cmsup\$(2)</p>``````

Result:
Check: fyk = 2 kN/cm2

Another interesting example is the following general design check routine:

``````#def check\$(x\$; x_max\$)
#if x\$ ≤ x_max\$
'The check is satisfied:'x\$'≤'x_max\$'✓
#else
'<p class="err">The check is not satisfied:'x\$'&gt;'x_max\$'✖</span></p>
#end if
#end def``````

Usage:

``````σ = 280MPa','σ_lim = 250MPa
check\$(σ; σ_lim)
τ = 120MPa','τ_lim = 0.58*σ_lim
check\$(τ; τ_lim)``````

Result:
σ = 280 MPa , σ lim = 250 MPa
The check is not satisfied: σ = 280 MPa > σlim = 250 MPa
τ = 120 MPa , τ  lim = 0.58·σ lim = 0.58·250 MPa = 145 MPa
The check is satisfied: τ = 120 MPaτ lim = 145 MPa

Modules, string variables and macros are parsed at preprocessor stage. Calcpad scans the code and rewrites the names with the actual contents. As a result, it generates “unwrapped” source code. If you check the respective checkbox you will see it in the output window, instead of the calculation results. It will also show automatically if any error occurs in preprocessing.

## 3. Angle units

°, deg – degrees;
– minutes = degrees/60;
– seconds = minutes/60;
rev – revolutions = 360 degrees.

You can convert between angle units and combine them with others to create composite ones like “rad/s” and “rev/min“. Angle units are accepted by all trigonometric functions. They override all other settings, made in the Calpad UI or by using “`#deg`“, “`#rad`” and “`#grad`” keywords in the code. Inverse trigonometric functions can optionally return the results with units, if you define the variable: “`ReturnAngleUnits = 1`“.

## 4. C/C++ style operators

Some relational operators in Calcpad use mathematical notation symbols like , , , , which are not available on computer keyboards. Although Calcpad provides buttons for insertion, this is still slower than typing. That is why, we introduced an alternative C/C++ like syntax that uses only keyboard symbols as follows:

The alternative symbols are replaced with the original ones while typing.

## 5. Calcs suppression with #noc

Previously, there were only two visibility control keywords regarding equations:
#equ – displays the equation and the result;
#val – displays only the final result without the equation.

Recently, we added one more option:
#noc – displays only the equation, without the result;

It not only hides the results, but also suppresses the calculations. This allows you to use Calcpad as a simple equation editor for faster typing of formulas. Unlike MathML and LaTeX, which use special syntax, Calcpad formats raw, simple text equations. At the end, you can export the entire content to MS Word or Libre Office with MathType/Math formatted equations.

## 6. Integration with Tanh-Sinh quadrature

Double exponential quadrature (Takahasi & Mori, 1974) is one of the most efficient numerical integration algorithms. In our tests, it demonstrated outstanding performance, compared to other popular methods: Integration-test-results.pdf. It is also adaptive and very accurate. The only problem is that the function must be continuous and smooth. However, due to its efficiency, we decided to include it as a second method, using the “`\$Integrate`” command.

The original algorithm has been improved further by Michashki & Mosig (2016) and Van Engelen (2022). In Calcpad, we made additional enhancements by precomputing and caching the abscissas and weights.

In all cases, you can continue using the old “\$Area” command, which is more general purpose. It is powered by the famous adaptive Gauss-Lobatto-Kronrod quadrature by Gander & Gautschi (2000).

## 7. Improved formatting

Sums, products and integrals are now formatted with the respective large symbol and limits on bottom and top. That gives the reports more professional feel and look. For example:

``S = \$Sum{(-1)^k/(2*k + 1) @ k = 0 : 1000}C(n; k) = \$Product{(i + n - k)/i @ i = 1 : k}I = \$Integral{x/(e^x - 1) @ x = 0.001 : 10}``

are formatted as follows:

## 8. Input fields defaults

In the previous versions, default values for input fields were specified by clicking the respective question marks in the code. When you saved the file, they were appended consequently at the end. This made uncomfortable using external editors like Notepad++ for writing Calcpad code. After introducing modules and macros, it was almost impossible to match the values with the corresponding input fields.

That is why, we changed the approach and made the default values to be added immediately after the question marks. You must enclose them with curly brackets, like: “`a = ? {2}`“. You can also override the values inside modules and macros by specifying the values as in the following example:

``````#include abc.cpd #{4; 5; 6}
d = a + b + c``````

If abc.cpd has the following contents:

``````a = ? {1}
b = ? {2}
c = ? {3}``````

The result will be:
a = 4
b = 5
c = 6
d = a + b + c = 4 + 5 + 6 = 15

## 9. Custom rounding

You can also override the global rounding settings for certain parts of your code, by using the “`#round` n” keyword. Replace n with the number of digits after the decimal point (from 0 to 15).

## 10. Pausing and step by step execution

If you have long calculations, you can pause the solution by pressing the “Pause/Break” key or “Ctrl + Alt + P” from the keyboard. Then, you can resume it with “F5“. Also, you can define breakpoints in your code by using the “`#pause`” and “`#input`” keywords. They correspond to two types of breakpoints:
`#pause` – stops the solution at the current line and displays the avalable results;
`#input` – stops at themcurrent line and renders an input form.

On the next `#input`, it will calculate and show the results before the current one and render the next input form. Sometimes, engineers need to see some intermediate results before entering more data. For example, when you design a RC beam, you calculate the required main reinforcement from bending ULS. After that, you select some actual reinforcement and use it further for SLS checks like deflections, crack widths, etc.

## 11. Performance and optimization

In the last few versions, we made important optimizations to improve Calcpad performance, as follows:

• UI responsiveness – it is not lagging anymore on large files;
• text processing and parsing – string splitting is replaced by the more effective memory spans wherever possible;
• function caching – now it works faster, but consumes much less memory;
• compiling – “constant folding” was introduced for faster execution;
• units – units arithmetic performance was enhanced by 20% – 30%.

## 5 reasons to choose math worksheets over Excel spreadsheets for structural calculations

Although there are plenty of software nowadays, engineers still have to develop custom calculations to solve specific problems. Some may still make them “by hand”, but most engineers use some kind of software for that.

The problem is that in most cases, we have to document our calculations and collect them in design reports or calculation notes. And it would be also very good, if we could give them professional look and feel. The other issue is to make the calculations reusable, in case we have similar projects in future.

Apparently, a peace of paper and pencil would not do the job anymore. Even a simple hand calculator (except for quick estimates). What we need is a suitable software for that purpose, but what?

To my experience, the most widely used solution for custom engineering calculations is Excel. At the same time, it is probably the most inappropriate for the job. But it is not actually an Excel fault. In general, it is a great software, very mature and high quality, but for tabular data. I myself usе it quite a lot for bills of materials and construction cost calculations. Spreadsheets are the best for simple formulas that have to be repeated on multiple rows/columns.

However, structural design algorithms are tree-like in nature, not tabular. Even if a tree data structure can fit in a table, this is not the optimal solution. It has some issues that cannot be neglected. During the years, I had to write a lot of structural design spreadsheets on Excel, as well as math worksheets on software like MathCAD or Calcpad. Over time, I found some important reasons to prefer the math approach over Excel spreadsheets:

## 1. Formulas are built by referencing cells.

It is slow to write and difficult to check and follow, especially for the longer ones. “M = FL + qL^2/2″ is much easier to type and more readable than “= D2A3 + C1A3^2/2″. It is like obfuscation of the actual formula. For longer calculations, it gets even harder. Imagine that you have 1000 lines (rows) of calculations and your formula uses some cells at row 1, 500, 200, 700, 900, etc. You will have to scroll up and down multiple times to find the required cells to be referenced. Of course, you can use named cells, to improve readability, but it requires additional work.

## 2. Formulas are hidden inside the cells.

That makes the printed report a real black box for the checking engineer. There is an option to make them visible, but without formatting they are difficult to read. The best option is to write them once again with MathType. Besides it is a double work, you can mistype something and there is no guarantee that the formulas you see, are the same inside the cells. Also, variable substitution in formulas is not available in this case. So, Excel spreadsheets are prone to errors.

## 3. You cannot apply units of measurement to values and use them in formulas.

You can write units of measurement in the next cell, but you cannot attach them to values and they will not take part in your calculations as it happens in math software. So, you must include the proper conversion factors manually in all of your formulas. Otherwise, you will simply get incorrect results.

## 4. It is difficult to define your own custom functions.

You must use VBA programming, which is not an option for a newbie. In math software, this is a piece of cake. You can simply type f(x) = 2*x^3 – … Also, it is very easy to plot a function graph directly from its equation. In Excel, you must calculate it in multiple cells and then, plot the respective range. If you have discontinuities or infinities and if you are not quite familiar with the function properties in advance, you can easily miss them and end up with incorrect graph.

## 5. It is almost impossible to branch the solution path and the report contents.

There are a lot of cases in structural design when the solution continues in completely different ways, depending on some intermediate checks. In systems like Calcpad, you can just add if-else if-else-end if condition blocks and put different formulas, text, and images in each one. In Excel, you can achieve limited results with intensive use of VBA.

So, even if you are an experienced spreadsheet developer and user, it is worth to give the math approach a try. Calcpad is perfect for this purpose, because it is lightweight, fast, easy to learn and free.

## How to calculate square/cubic root by hand?

Today, technologies make everything easy and fast. And sometimes, when we see the ancient marvels of engineering, we wonder: “How they designed this without computers”? And the answer is: “with knowledge, inspiration and hard work”.

However, being almost irreplaceable, computers makes us more and more adjective. Our brains become lazy and we are losing some knowledge and skills. Such is the manual calculation of root. It is so easy and elegant, so I will share it here for the sake of the good old math science.

## Square root

Let’s calculate x = √a first. By squaring both sides and moving a to the left side, we obtain the following equation:

x2a = 0

This is an equation of type f(x) = 0 and we can solve it numerically, using the Newton’s method, as follows:

1. Start with initial guess: x0. It is good to be close to the solution, so you can use the nearest exact root.
2. Calculate the next approximation of the root by the equation: x1 = x0f(x0)/f′(x0)
3. Continue the iterations by using the previous result to estimate a new one:
` `xn+1 = xnf(xn)/f′(xn).
4. Stop when you reach the desired precision.

In the above formulas, f(x) is the function and f′(x) is the first derivative. In our case:
f(x) = x2a
f′(x) = 2x

If we substitute the above equations into the Newton’s iterative formula, we get:
xn+1 = xn − (xn2a)/2xn
xn+1 = (xn + a/xn)/2

Now, lets use it to calculate √a = ?, for a = 5. We will start with an initial guess of x0 =√4 = 2.

1st iteration: x1 = (x0 + a/x0)/2 = (2 + 5/2)/2 = 2.25
2nd iteration: x2 = (x1 + a/x1)/2 = (2.25 + 5/2.25)/2 = 2.2361
3rd iteration: x3 = (x2 + a/x2)/2= (2.2361 + 5/2.2361)/2 = 2.2361

After the second iteration, we have 5 correct significant digits, which is quite sufficient for most cases.

## Cubic root

We will use the same approach, as follows:
f(x) = x3a
f′(x) = 3x2
xn+1 = xn − (xn3a)/3xn2
xn+1 = (2xn + a/xn2)/3

Let’s see how it works for 3√7. Our initial guess is: x0 = 3√8 = 2.

1st iteration: x1 = (2x0 + a/x02)/3 = (2·2 + 7/22)/3 = 1.9167
2nd iteration: x2 = (2x1 + a/x12)/3 = (2·1.9167 + 7/1.91672)/3 = 1.9129
3rd iteration: x3 = (2x2 + a/x22)/3 = (2·1.9129 + 7/1.91292)/3 = 1.9129

Again, we managed to find it after two iterations. We can also use Calcpad to create a tiny program for that:

It is interesting that after just 4 iterations, the Newton’s method finds 16 correct significant digits, which is the limit precision for most computers. We also calculated the root by three different methods which gived the same result.

## What’s new in Calcpad, version 5.6.3?

Since version 5.6 was released, we made some little, but nice and useful improvements. They were driven mostly by our small and dedicated GitHub community. Some of the most important changes are listed below:

### 1. Autorun mode

When enabled, the results are refreshed automatically each time you edit the code and leave the current line. You can switch it on and off by the checkbox over the output window. When the “Autorun” is off, the program displays “Quick help” on the right as before. You can run calculations by pressing the respective button. However, new shortcuts were added to make that easier. Press F5 to run the calculations and F4 to compile to input form. Ctrl + Enter calculates the results and scrolls the output to match the current line.

### 2. Scroll enhancements

After recalculation, the scroll position of the output window is preserved and not moved to the beginning as before. However, if you get lost, you can double click on a line in the source code and the output will scroll to match. Alternatively, you can press Ctrl + Enter.

### 3. New symbols

The Symbol toolbar was extended with the full set of Greek letters. You can also change the case by pressing the button on the right.

There is also a shortcut to convert Latin to Greek letters. You can use it to insert Greek letters even faster. Just type the corresponding Latin letter and press Ctrl + G. The following conversion table is used:

Some special symbols were also added: ” ø “, “Ø”, ” ° “, “∡” and primes: ” ′ “, ” ″ “, ” ‴ “, ” ⁗ “. You can now use them in variable and function names. However, you cannot start a name with a number of prime symbol. Although primes look similar to quotes, they are different characters. So, there is no danger for the program to confuse them and start a comment.

### 4. Bugs and errors

There was a problem with recognition of inches (in) units in Complex mode that is now fixed. The “i” character was parsed as the imaginary symbol first and “n” dropped out as invalid unit.

Rewriting of non-metric (Imperial and USCS) units was also improved. For example, in the previous versions, lbf/in^2 was rewritten as kPa by default. Now, ksi is used instead. It is the same with other mechanical units.

Some rare bugs were luckily identified and fixed as well. For example, if you attempted plotting with identical limits, the program fell in infinite loop.

Debugging of Calcpad programs was made easier. Links to the respective lines in the source code were added to the error messages in the output. Also nice background highlighting was implemented.

Notepad++ is one of the most popular and free text/code editors. It supports many different programming and scripting languages. You can also edit HTML, XML, CSS, JSON and other types of files. Currently, its text editing capabilities significantly exceeds those of the Calcpad own editor. Some features that worth mentioning are:

• multiple file tabs;
• predefined and custom styles;
• column editor;
• auto completion;
• bracket matching;
• line operations/sorting;
• tab/white space management;
• content folding, etc.

And and last but not least – syntax highlighting. Unfortunately, the Calcpad language is not natively supported by Notepad++. However, there is an option to add a custom language with user defined syntax (UDL). Then, you can export it as an xml file and import it everywhere else.

The good news is that we have already created the required language definition files, so you can use them to write Calcpad programs with Notepad++.

To use Notepad++ as Calcpad code editor you need to the following:

3. Unzip the archive in the preferred folder. You will find two files inside:
– Notepad++.xml – for light mode and
– Notepad++dark.xml – for dark mode.
5. Click the “Import…” button in the dialog that will appear.
6. Browse for the one of the xml files above and click “Open“.

From now on, each time you open a .cpd file, it will be automatically displayed with the Calcpad syntax highlighting. For best results with the dark theme, select the appropriate style from the “Style Configurator” dialog.

You can find it in the “Settings” menu. For example, a good option is the Twillight theme. You can also set a nice monotype font, such as “Noto Mono“, “Liberation Mono“, “Source Code Pro“, “Consolas“. etc. You can also make the whole interface t dark, by checking “Enable dark mode” in the “Preferences” dialog from the same menu.

You can also make possible to run the Calcpad programs directly from inside the Notepad++ editor. For that purpose, you will have to install the Franco Stellari’s RunMe plugin for Notepad++. You can find detailed instructions in this file: install.txt. It will add a small button at the end of the Notepad++ toolbar . Click the button or press Shift+F5 to run the current program.

.

## How to plot the Julia Set

In the previous post we discussed the Mandelbrot set. There is another beautiful algebraic fractal, called Julia set 𝓙. It is named after the French mathematician Gaston Julia, who has discovered. It is obtained by the same formula as the Mandelbrot Set:

zn+1 = zn2 + c

The difference is that z0 = x + yi is the current point (x, y) in the complex plane and c is a fixed complex number, while for the Mandelbrot set z0 = 0 and c = x + yi.

Both sets are closely related. There are a lot of different Julia sets, corresponding to different values of c. Their shapes and properties depend on where the point is located, related to the Mandelbrot set 𝓜. Inside 𝓜, the corresponding Julia set 𝓙 is connected. Outside 𝓜, it falls apart into a set of isolated points, called Fatou set 𝓕. Further away, it turns into Fatou dust – infinite number of points with zero area, similarly to the Cantor dust.

For simple plotting, you can use the following code:

```"Julia set
'Define the function
JuliaSet(z; c) = \$Repeat{z = z^2 + c @ i = 1 : n}
#hide
'Set the plotting parameters
PlotStep = 1','PlotWidth = 600','PlotHeight = 400
#show
'Plot for'c = -0.4 + 0.59i'and'n = 200'iterations
\$Map{abs(JuliaSet(x + 1i*y; c))) @ x = -1.5 : 1.5 & y = -1 : 1}
```

When you run it inside Calcpad, you will get the following result:

To obtain better and colorful images, we will use different approach this time. We will modify the JuliaSet function in a way to return directly the number of iterations for which the result overpasses 2. The source code is provided bellow:

```"Julia set
J(z; c; n) = abs(\$Repeat{z = if(abs(z) < 2; z^2 + c; 4) @ i = 1 : n})
JuliaSet(z; c) = \$Find{J(z; c; x) - 2 @ x = 1 : n}
#hide
PlotStep = 1','PlotWidth = 600','PlotHeight = 400','Precision = 10^-6
#show
'Plot for'c = -0.4 + 0.6i'and'n = 100'iterations
\$Map{JuliaSet(x + 1i*y; c)) @ x = -1.5 : 1.5 & y = -1 : 1}
```

It runs a little bit slower, but the result is worth the wait:

You can download Calcpad for free and experiment by yourself with different values of c. For example, the lightning image bellow is a special case of Julia set for c = i and is called “dendrite” fractal.

## Plotting the Mandelbrot Set

Calcpad has an interesting and undocumented feature, that I am going to reveal in this post. You can use it to quickly plot the Mandelbrot set. This is a set of complex numbers c, for which the iterative equation zn+1 = zn2 + c does not go to infinity. The most beautiful part is that it shows fractal behavior at its boundaries. If you plot it with the appropriate colors, you can create stunning images and animations. I really like those on the Maths Town’s YouTube channel:

Usually, such videos are created by complicated software, specially developed for that purpose. However, you can also plot the Mandelbrot set with Calcpad and a few lines of code. There are two ways to do that. The first one is to use the Calcpad capabilities for complex arithmetic, as follows:

```"Mandelbrot set
'Define the function
MandelbrotSet(z; c) = \$Repeat{z = z^2 + c @ i = 1 : n}
#hide
'Set the plotting parameters
PlotStep = 1','PlotWidth = 500','PlotHeight = 500
#show
'Plot for'n = 50'iterations
\$Map{abs(MandelbrotSet(0; x + 1i*y)) @ x = -1.5 : 0.5 & y = -1 : 1}
```

Open Calcpad, switch to “complex”mode, paste the above code inside, and run the calculations. You will get the following result:

It is not too bad, but still not very impressive. At least, the exterior is missing any colors. Then, how we can do the coloring? Theoretically, it is satisfied for all numbers in the set that |zn| < 2 for any n. Outside, the process diverges and the modulus goes to infinity, as n increases. You can use that in the following way:
1. For each point, calculate the number of iterations n, where |zn| gets larger than 2.
2. Map these values to colors.
3. Plot the colors at the respective coordinates.

If you play with the number of iterations n in the above example, you will see that the plot changes. The shaded area is where |zn| converges for the respective n.

If you combine the above plots and apply some colors, you can get really beautiful artistic image like the one bellow:

The only problem with this method is that it is slow and elaborate. That is why, I included a special function Mandelbrot(x; y) in Calcpad. It calculates the number of iterations, needed for the modulus to overpass 2. You can plot this function directly with the preferred color scale:

To do that, I used even shorter piece of code:

```"Mandelbrot set
#hide
'Set the plotting parameters
PlotStep = 1','PlotWidth = 500','PlotHeight = 500
#show
'Plot for'n = 50'iterations
\$Map{Mandelbrot(x; y) @ x = -1.5 : 0.5 & y = -1 : 1}
'Plot for'n = 70'iterations
\$Map{Mandelbrot(x; y) @ x = -0.59 : -0.58 & y = 0.55 : 0.56}
```

Here you can find the C# code of the Mandelbrot function that works inside Calcpad:

```        private static readonly double Log2Inv = 1 / Math.Log(2);
//Calculates the number of iterations for which z = z^2 + c satisfies |z| > 2
//Then transforms it to a smooth log scale and returns the result
protected static double MandelbrotSet(double x, double y)
{
//Checks if the point is inside the set and directly returns NaN
if (x > -1.25 && x < 0.375)
{
if (x < -0.75)
{
if (y > -0.25 && y < 0.25)
{
double x1 = x + 1,
y2 = y * y,
x2 = x1 * x1;
if (x2 + y2 <= 0.0625)
return double.NaN;
}
}
else if (y > -0.65 && y < 0.65)
{
double x1 = x - 0.25,
y2 = y * y,
q = x1 * x1 + y2;
if (q * (q + x1) <= 0.25 * y2)
return double.NaN;
}
}
//For all other points performs detailed calculations
double re = x, im = y;
for (int i = 1; i <= 1000; ++i)
{
double reSq = re * re,
imSq = im * im,
sumSq = reSq + imSq;
//To avoid the sqrt function, the check |z| > 2 is replaced by z^2 > 4
if (sumSq > 4)
{
var logZn = Math.Log(sumSq) / 2;
var nu = Math.Log(logZn * Log2Inv) * Log2Inv;
return (1.01 - Math.Pow(i - nu, 0.001)) * 1000;
}
//Calculates z = z^2 + c
im = 2 * re * im + y;
re = reSq - imSq + x;
}
return double.NaN;
}
}
```

The complete source code is available on GitHub.