Arduino Dice

In the previous post I outlined the hardware of the multifunction shield. The plan now is to build a few simple projects based around it, so lets start with something really simple. We have a button, and we have a screen, so lets generate some random numbers and make an electronic dice.

make button1 digital input A1

make segLatch digital output 4

make segClock digital output 7

make segData digital output 8

make segVal number

when shiftOut

.repeat 8

..if segVal > 127

...set segData to on

...set segVal to segVal-128

..else

...set segData to off

..

..set segVal to segVal*2

..

..set segClock to on

..set segClock to off

make digit number

make digitVal number

make decimalP boolean

when showDigit

.set segVal to 255

.if digitVal=0

..set segVal to 192

.if digitVal=1

..set segVal to 249

.if digitVal=2

..set segVal to 164

.if digitVal=3

..set segVal to 176

.if digitVal=4

..set segVal to 153

.if digitVal=5

..set segVal to 146

.if digitVal=6

..set segVal to 130

.if digitVal=7

..set segVal to 248

.if digitVal=8

..set segVal to 128

.if digitVal=9

..set segVal to 144

.

.if decimalP

..change segVal by -128

.

.set segLatch to off

.broadcast shiftOut and wait

.if digit=1

..set segVal to 241

.if digit=2

..set segVal to 242

.if digit=3

..set segVal to 244

.if digit=4

..set segVal to 248

.broadcast shiftOut and wait

.set segLatch to on

.wait 2 millisecs

Everything above is standard code which sets up the 7 segment display on the multifunction shield. Show digit allows us to display a value on any of the segments, with an optional decimal point.

Normally we'd also include code to display the number in the string message on the digits of the display but as we just want to display a single digit, we can simplify (and make it look more fun) by displaying the same value on all four digits:

when start

.set decimalP to no

.forever

..repeat 4 using digit

...broadcast showDigit and wait

Now to generate a random number every time we press the button:

when start

.forever

..wait until not button1

..set digitVal to pick random 1 to 6

..wait 0.1 secs

..

..wait until button1

..wait 0.1 secs

We wait until the button is pressed (logic low), and set the message to the value of the number.

After the button is pressed and we've updated the digitVal, we wait 0.1 seconds then wait for the button to go high again. This is to "debounce" the switch - when you press a switch it can bounce a little, so flips high and low for a fraction of a second before being firmly pressed. Adding a little delay stops this happening.

Now try it:

2,2,6,3,5,3,1...

Never gamble with a computer! There's a serious problem. Every time you run the code you get exactly the same results. Well its a computer after all - programming depends on computers doing the same thing over and over again. When you ask a computer to generate a random number, it runs some code to "make something up", but every time that code is run it does the same thing, so we always end up with the same sequence!

In fact its really hard to generate random numbers with a computer. Why do you think the lottery uses balls with numbers written on them, rather than something "smarter"? There are actually companies that will sell you random numbers (they let you have them free for non-commercial use!).

We need to add some genuinely random element to the system - like a user. Here's an alternative version of the code:

when start

.forever if not button1

..set digitVal to pick random 1 to 6

..wait 0.01 secs

We pick a new random number ever 100th of a second until the user lets go of the button. As the numbers are changing to fast for us to see, and we can't accurately move/remove our finger from the button, this  is pretty random.

But how could we test it?

If you were just interested in a little light hardware hacking, you should probably stop reading now... I'm about to go back into a whole load of maths!!! Sorry... bad habit.

Well we could start by throwing a sequence of numbers and seeing what it looks like:

3,5,6,3,3,1,1,5,6,6,1,6,5,6,1,1,2,3,1,3,5,3,3,4

I generated 24 numbers, which means we'd expect to see 4 of each number:

1:6

2:1

3:7

4:1

5:4

6:5

Not even close! I got hardly any 2's or 4's. What's the chance of rolling a dice 24 times and only getting one 4 on the very last go? Well the chance of not getting a 4 on the first go is 5/6, on two rolls it's 5/6*5/6 and so on, which means that the chances of getting 23 non-4's followed by a 4 would be (5/6)^23 * 1/6 or about 0.25% - not very likely.

But not having many fours is just one of the many thing we might notice. There are 6 values on the dice. Any go them only coming up only once would be equally suspicious. Or not coming up at all. or coming up as the first and last number, or coming up in a massive streak in the middle - there are lots of "interesting" patterns we could spot. In fact perhaps the most "interesting" thing that could happen is for nothing interesting to happen. If I didn't see at least something suspicious I'd be very suspicious!

This is sometimes called the Texas sharp shooter fallacy - after a guy who fires a load of bullets into the side of a barn, then walks over and paints a target around the largest cluster. We need to say what we're looking for prior to collecting the data. Here we've seen something interesting in the data, and calculated the chance of it happening - except of course the real chance of it happening is 1. IT ALREADY HAPPENED. We started with a sequence with only one 4. Thats what made us look for only one 4, so of course we found only one 4!!!!

So you do we tell if that sequence is dodgy? To do it properly we need to start with our assumption (called the null hypothesis), and then look for evidence that this is not true! In this case our null hypothesis is that the dice is fair, and for any number is equally likely on only given throw. Then, we throw the dice, and then we do a Chi Squared test.

OK - that looks scary. But lets breakout down. Ei is the number of times we expect a number to appear - we threw 24 times, and have 6 sides, so Ei=4. Oi is the number is the number of times it actuallyy happened (Observed), so for each side we calculate (O-4)/4 and then add them up.  Heres the calculation for the previous run:

1:6: (6-4)^2/4=1

2:1: (1-4)^2/4=2.25

3:7: (7-4)^2/4=2.25

4:1: (1-4)^2/4=2.25

5:4: (4-4)^2/4=0

6:5: (5-4)^2/4=0.5

Adding these gives us chi squared value of 8.25. We then compare this to a value from a table. The only catch is that there are several versions of the table depending on the "degree's of freedom". Degrees of freedom is just a fancy way of saying how many different things can change and in this case out dice has 6 sides, so we have six values that can change... BUT the last side doesn't count! Why? Because if I tell you I threw 24 times, and I tell you how many times I got the values 1-5 you already know the final answer, so the degree's of freedom is actually only 5.

Degrees of freedom	Probability less than the critical value
	0.90	0.95	0.975	0.99	0.999
5	9.236	11.070	12.833	15.086	20.515

Here's the table we need to look things up in. If the number we get is bigger than the one in the table, then our "null hypothisis"is unlikely. If we get value of greater than 11 then there's only a 5% chance of getting that result if the dice is fair, and we "reject the null hypothosis".

We got a value less that 9.2, so its fair right? No!!!!! We haven't shown its fair - just that we haven't got enough to prove that its NOT fair. Even worse - even if it is fair we'd still think it was unfair 10% of the time!

So lets try and collect more evidence.

make counts list of numbers

make trials number

make rolled number

make tmp number

make expected number

make chi number

when rollFairDice

.delete all of counts

.repeat 6

..add 0 to counts

.

.repeat trials

..set rolled to pick random 1 to 6

..set tmp to item rolled of counts

..change tmp by 1

..replace item rolled of counts with tmp

when start

.

.set trials to 10*6

.say join "Trials:"[trials]

.

.broadcast rollFairDice and wait

.

.set expected to trials/6

.say join "Expected:"[expected]

.set chi to 0

.repeat 6 using rolled

..say join join [rolled] ":" [item rolled of counts]

..set tmp to item rolled of counts

..set tmp to tmp-expected

..set tmp to tmp * tmp

..set tmp to tmp/expected

..change chi by tmp

.

.say join "chi:"[chi]

This code "rolls the dice" 60 times (so we'd expect to get an average value of 10 for each face), and on my Mac produces:

1:7

2:14

3:13

4:9

5:9

6:8

chi:4

600 trials

1:102

2:89

3:121

4:101

5:90

6:97

chi:6.76

6000 trials

1:980

2:993

3:1030

4:1009

5:1002

6:986

chi:1.63

We can see that the chi squared value is changing but not in any obvious way - different values of trials can produce some bigger values. Remember that just by chance, we'll get a value of 15 or more 1 time in 100, no matter how many trials we have!

What lots of trials will to is make that number go super big by accumulating evidence, if there is anything fishy going on. Here we're looking for evidence of impropriety when in fact there is none, which is pretty dull... so lets start cheating:

when rollBiasDice

.delete all of counts

.repeat 6

..add 0 to counts

.

.repeat trials

..set rolled to pick random 1 to 7

..if rolled=7

...set rolled to 6

..set tmp to item rolled of counts

..change tmp by 1

..replace item rolled of counts with tmp

Now 6 will come up twice as often as any other. Lets roll 24 numbers like we did before:

2,5,3,6,2,4,4,3,2,4,3,6,6,5,6,1,4,2,3,3,4,1,6,3

Expected:4

1:2

2:4

3:6

4:5

5:2

6:5

chi:3.5

It really doesn't look that different! Chi squared of 3.5 means the test pulled up nothing suspicious, so we're all good right? NO! we just haven't got enough evidence yet. We KNOW we're cheating, and 6 has come up more times than we'd expect, but the nature of probability is that its very likely that this could happen by chance (in fact 3 has come up more than 6, and that is purely by chance).

Lets do 120 rolls:

Expected:20

1:20

2:17

3:19

4:18

5:15

6:31

chi:8

Here we're starting to see 6 pull clearly ahead, and the chi squared  is reflecting that, but statistically we're not there yet. If we're running a casino, then we're looking out for lots of 6's, but that's not what we're testing here - a high or low score of any number would be sufficient to reject the null hypothisis that the dice is fair, so there are many ways we could get a result this "interesting".  Maybe we should have used the null hypothosis "6 appears with a probability of 1/6" but that would have mean different numbers for the test, which may produce a higher chi squared. Remember the sharp-shooter: we're not allowed to change the test after we've seen the results!

But look what happens when we put the number of trials up to 240:

Expected:40

1:38

2:31

3:33

4:38

5:39

6:61

chi:14.5

Chi squared keeps going up because now we really are collecting more evidence - but looking at the table, this (or something "equally interesting") could still happen more than 1 time in 100.

Going to 480 trials finally kicks us over the edge:

Expected:80

1:68

2:61

3:81

4:74

5:75

6:121

chi:28.1

I think the take away here is that you can't prove the null hypothosis - only collect evidence to reject it. When we ran our tests with a fair dice, nothing much happened. Increasing the trials didn't change anything - if all is fair, then chi squared will be bigger or smaller at random - we keep testing, but we never find more or less evidence, because there's none to find. We can't produce evidence that the dice is unfair, because there's no evidence to find! However when we did have some thing dogey going on, increasing trials gave us more evidence.

There's a guy who takes this pretty seriously and built a machine to roll dice, and then read the values of with a camera, and image recognition system!

This really was supposed to be a fun/simple hardware build, but I guess it got side tracked, but I think it was worth it. You really should simulate a bias dice (make it more or less bias), and see how long it takes to detect. Running stuff like this in code helps you get a much better feel for the maths.

Also while Chi Squared is A level maths, the actually implementation is pretty easy. Implementing this sort of code is something you could do at KS4, and it would make a lot of sense if tied to some kind of practical work.

Arduino Multifunction shield

In preparation for a robot workshop I ordered a bunch of Arduino Motor Controller Shields. These are super cheap from china, and can drive up to 4 DC motors. They also break out a couple of pins on servo headers, so they're ideal for getting kids to actually assemble and program a robot.

A few weeks later a brown box arrived with the expected customs declaration "other electronic components". These turn up quite regularly here at Sniff Manor, and with the workshop still a month or so away I put the box to one side "for later inspection". As a bunch of "more interesting" stuff turned up, the box sat for a few days without attracting much attention, before I figured I should put it away in the big box of parts I'll need for the workshop. However before putting it away I thought I best check that it was in fact the motor controllers - I do sometimes order stuff and then forget I ordered them...

In this case it was just as well as while the box did contain Arduino shields, they weren't the one I expected. Someone working on a packing line several thousand miles away, had sent me the correct number of the wrong part. They were superficially similar but completely not what I wanted.

A quick search revealed that I'd been sent a bunch of "Multifunction" shields. I then quickly ordered a new set of motor shields (should be here in time!), before deciding there really was no point in sending back the multifunction shields - for what the whole batch cost it wasn't worth the hassle of trying to get them replaced with the right thing. Besides - a new component to play with!!!

The Multifunction shield is a but of a mash up (as its name suggests), which I guess is aimed at beginners wanting quick and easy way to hook up a load of components and start programming them. Unfortunately like many Arduino shields, documentation is scarce, but some searching did turn up an circuit diagram, and some sample code - all written in Chinese. Working through these, and a bit of poking with a multi-meter, I was I was able to get the board up and running.

make led1 digital output 13

make led2 digital output 12

make led3 digital output 11

make led4 digital output 10

make button1 digital input A1

make button2 digital input A2

make button3 digital input A3

when start

.set led1 to on

.set led2 to on

.set led3 to on

.set led4 to on

.

.forever if not button1

..set led1 to off

..wait 0.2 secs

..set led1 to on

..

..set led2 to off

..wait 0.2 secs

..set led2 to on

..

..set led3 to off

..wait 0.2 secs

..set led3 to on

..

..set led4 to off

..wait 0.2 secs

..set led4 to on

The first and simplest things to tackle are the 3 switches and  4 LEDS. Which are attached to the pins of the Arduino, as outlined in this code. The switches have pull-up resistors attached (which you can disable by removing J2, though that still leaves the pins connected together by 10K resistors, and in any case the pins aren't broken out elsewhere, so really not much use). However they are pull-UP's, which mean that pressing the switch generates logic low. Similarly the LED's are also active low - not ideal in a board that's supposed to appeal to newbies: set led to on, turns it off, and set led to off turns it on! There's no reason to design a board this way, and while many "real" circuits are active low its enough to make me hesitate in using this with kids. 

make buzzer digital output 3

when start

.forever

..if not button3

...set buzzer to off

..else

...set buzzer to on

The next component we'll look at is the buzzer. Note that this is a buzzer, which makes a very loud, fixed pitch beep, rather than a piezo transducer found on the pibrella. That means its easy to code - we just run it on and off using pin 3 (though again its active low).

make pot analog input A0

when start

.forever

..set message to [ pot *10 ]

..wait 0.1 secs

The blue rectangular block is a 10k pot, set up as a potential divider on A0. While pots can make great user input devices, this one uses a tiny screw, which requires many turns to change it. It's handy for demonstrating the idea, but something a bit bigger, and perhaps with less turns would have fitted the purpose of demo'ing analog input a bit better.

make segLatch digital output 4

make segClock digital output 7

make segData digital output 8

Perhaps the most useful feature of the board is the set of 4x7 segment displays. These are attached to three pins: clock, data and latch. These are the sort of thing we could write a "device" for, but its actually not that hard to program them directly in Sniff.

make segVal number

when shiftOut

.repeat 8

..if segVal > 127

...set segData to on

...set segVal to segVal-128

..else

...set segData to off

..

..set segVal to segVal*2

..

..set segClock to on

..set segClock to off

This is the low level code which loads a byte into the 7 segment display registers. To display a character we first load in a byte showing which segments should be illuminated, and then a second byte turning on (usually) one of the digits. Only one of the digits is actually lit up at a time, but we can keep cycling and it will look fine.

make digit number

make digitVal number

make decimalP boolean

when showDigit

.set segVal to 255

.if digitVal=0

..set segVal to 192

.if digitVal=1

..set segVal to 249

.if digitVal=2

..set segVal to 164

.if digitVal=3

..set segVal to 176

.if digitVal=4

..set segVal to 153

.if digitVal=5

..set segVal to 146

.if digitVal=6

..set segVal to 130

.if digitVal=7

..set segVal to 248

.if digitVal=8

..set segVal to 128

.if digitVal=9

..set segVal to 144

.

.if decimalP

..change segVal by -128

.

.set segLatch to off

.broadcast shiftOut and wait

.if digit=1

..set segVal to 241

.if digit=2

..set segVal to 242

.if digit=3

..set segVal to 244

.if digit=4

..set segVal to 248

.broadcast shiftOut and wait

.set segLatch to on

.wait 2 millisecs

This code shifts out the two bytes representing  a number (from 0-9) and displays it on the selected digit, optionally setting the decimal point.

make message string

when start

.make index number

.make doneDP boolean

.forever

..set index to 1

..set doneDP to no

..repeat 4 using digit

...set digitVal to value of letter index of message

...set decimalP to no

...if index=length of message and not doneDP

....set decimalP to yes

...change index by 1

...if letter index of message ="."

....set decimalP to yes

....set doneDP to yes

....change index by 1

...broadcast showDigit and wait

Finally here's the code that goes through the string message, and sends 1 digit of it at a time to the display. The tricky part here is handling the decimal point. If we're displaying a number, and the next character is a decimal point then we display the point with this digit, and skip over the next character. When we get to the end of the string we print a point if we haven't already, and pad with zero's.

In the C examples that I found for  the board the code has to keep calling display in the main loop to keep the lcd updated - remember only one set of segements is lit at once, so we need to constantly keep refreshing the display.  Fortunatly Sniff can handle this more easily by just running this script when the program starts. It keeps going forever, updating the display in the background. Earlier we wrote a script which reads a value from the pot and assigns it to the string message. Now this is constantly displayed on the LED display.

In fact because the Sniff scripts that make up a program can all run at the same time, we can run the code for the pot, the display, the leds, and the buzzer all at the same time, in a single demo, rather than writing a demo for each feature!

In addition to the built in devices on the boardthere are breakout pins and sockets intended for various purposes.

make servo1 digital output 5

make servo2 digital output 6

make servo3 digital output 9

make servo4 digital output A5

There are four 0v/5v/Signal servo headers. Of course to use them with servo's you would make them digital outputs, but they can equally be used as inputs (and/or analog if the pin supports it).

The block of 7 header in the top left has 0/5v on the first two pins, and serial (Arduino pins 0 and 1) on two of the others.

make irSensor receiveIR device 2

#make lm35 analog input A4

make thermometer ds18 device A4

The block of 6 pins across the middle is designed of attaching an IR receiver on pin 2 in the first three holes, and a temperature sensor on A4 in the second. From left to right the pin out is D2,0V,5V (for theIR receiver) then 0V,A4,5V for the temperature sensor. The temperature sensor can be an LM35, or a ds1820. In the case of the DS1820 J1 provides the required 10K pull-up on the data line - it can/should be removed it you're using A4 for something else. Both a TSOP4838 IR receiver, and a DS1820 should just plug straight in. Though its not obvious there are little white outlines on the board - match the shape, to know which way around.

So that's the technical spec's of the board. The buttons, display, and servo header pins are all pretty useful,  as together these have the potential to work well for a range of experiments, the negative logic is annoying in a beginners board.  It's also a shame that i2c isn't broken out. I know that's advanced for the target audience, but this would give access a bunch of sensors, which could simply be read, and displayed on the 7 segment display - instant win! A4 and A5 are broken out and unused, but in different parts of the board, so it would be messy (you could install a breakout/io shield underneath just to get i2c, but really?).

I could definitely see then being useful in a workshop. Give me a few days, and I'll post some examples of using the board to do something fun!

Beyond Newton Rhapson

Normally I blog about fun programming exercises you can do using Sniff. Often this involves robots, flashing lights, and physics experiments. They're usually appropriate for kids somewhere in the age range 10-16. Recently however we've got a little off track. It started innocently enough with a simple Monte Carlo Simulation, but then moved on to root finding. It got a bit out of hand with Taylor Approximations, and then Newton Rhapson. If you're looking for a fun exercise for you KS4 physics class check the older posts on the right!

In the meantime here's one last maths heavy post (promise!)

Newton's Method for root finding (or strictly the Secant Method, as we're finding the derivative analytically) works pretty good. It fits a straight line to a function, then solves for the straight line. Then we repeat again using the straight line solution as a starting point to fit a new straight line.

We're about to leave anything that's relevant to anything that's taught in schools - so lets abandon any pretence this isn't about calculus. I'll also abandon any claim that I actually know this stuff - my maths is pretty good, but I don't have a maths degree, so there may be some inaccuracies in here! This is stuff I picked up when I was fact checking the previous post!

Newton Rahpson takes the derivative of a function find its gradient at a point and fit a straight line, but what if we could do better? The function isn't a straight line - its a curve, which means its gradient is changing, so we get the wrong answer. But what if we could include information about not just the gradient, but how the gradient was changing? Well if we plotted a graph of the gradient, what would be it's gradient? In other words the second derivative!

For a straight line it has a constant gradient so f'(x)=m which is a constant. f''(x) is how the gradient changes, and (for a straight line) it doesn't so f''(x)=0. But if we had a quadratic (say x*x), then the gradient does change - its a curve! In fact in this case f'(x)=2x. That of course is a straight line with gradient 2, so f''(x)=2. Newton assumes that f''(x)=0 which is why it only approximates the curve rather than getting it right.

It turns out Newton's method is only the most basic approach to this problem and there are a whole bunch of more advanced solutions which incorporate higher order derivatives, and fit ever more complex curves. After Newton, the next simplest method is Halley's method. It's just like Newton Rhapson except that we use:

That's a bit more complex, but most of it is just an equation that we need to type in, in place of the Newton equation. As you can see it has f''(x) in the bottom right, so its accounting for the curvature.

We approximate f'(x) as:

   f'(x) = (f(x+dx)-f(x))/dx

so we can approximate f''(x) as:

   f''(x) = (f'(x+dx)-f'(x))/dx

However using our equations so far this would approximate f'(x+dx) using f(x+dx) and f(x+dx+dx), which is a bit far away from where we're interested in. However we could just as easily have used:

   f'(x) = (f(x)-f(x-dx))/dx

(this is called backward differencing, rather than forward differencing), so we can use these two approximations - one backwards and one forwards to get two derivatives close to x, and use those to find f''(x) so:

f''(x)=((f(x+dx)-f(x))/dx - (f(x)-f(x-dx))/dx )/dx

It's the forward difference, minus the backward difference, divided by the distance between them. Rearranging this gives:

f''(x)=(f(x+dx)+fx(x-dx)-2f(x))/(dx^2)

(useful info here)

Also while we're at it we can get a better approximation for f'(x) by combining both the forward and backward differencing equations to get:

   f'(x) = (f(x+dx)-f(x-dx))/2dx

In fact we can drop this into Newtons method for slightly better results.

Now we can find f''(x), and a better approximation to f'(x) we can just drop these into the previous code, and off we go:

when start

.set x to 45

.set dx to 0.005

.repeat 10

..broadcast calcF and wait

..set fx to f

..say join join [x] ":" [fx]

..

..change x by dx

..broadcast calcF and wait

..set fpdx to f

..change x by -dx

..change x by -dx

..broadcast calcF and wait

..set fmdx to f

..change x by dx

..

..set m to (fpdx-fmdx)/(dx+dx)

..set d2 to (fpdx+fmdx-(fx+fx))/(dx*dx)

..

..set x to x-((2*fx*m)/(2*m*m-fx*d2))

Running that on our random cubic, prints out:

45:95129

29.7889:28177.9

6.22419:311.385

3.00481:41.1829

1.425:4.52992

0.869719:0.30097

0.802179:0.00100017

0.801938:-1.19209e-07

0.801938:-1.19209e-07

Compare that to the results from Newton:

45:95129

29.5867:27619.5

19.4132:8049.59

12.7117:2363.52

8.29118:698.161

5.37315:206.495

3.45087:60.4607

2.20288:17.1924

1.43379:4.62523

1.00966:1.05843

0.835976:0.145962

0.803114:0.00487411

0.801939:7.15256e-06

0.801938:2.38419e-07

0.801938:-1.19209e-07

Including the information about the curvature get us to the first answer a lot quicker! Using -1 and -2 as starting points to find the other roots produces:

-1:1

-0.599916:0.103806

-0.554994:8.11815e-05

-0.554958:1.19209e-07

-0.554958:-5.96046e-08

-0.554958:-5.96046e-08

-2:1

-2.23067:0.0828772

-2.24698:1.66893e-05

-2.24698:-9.53674e-07

-2.24698:7.15256e-07

And for Newton:

-1:1

-0.500083:-0.124813

-0.555568:0.00140119

-0.554958:-5.96046e-08

-0.554958:-5.96046e-08

-2:1

-2.33341:-0.481945

-2.2531:-0.0317504

-2.247:-0.00013113

-2.24698:7.15256e-07

-2.24698:-9.53674e-07

Both of these converge pretty quick, so the results are a little inconclusive: for -2 Halley converges 1 step faster. For the root near -1 it looks like Newton gets their first... 4 steps rather than 5, but look more closely:  on the first three iterations Halley is closer, and on the fourth is actually printing the right answer, but the error is 10e-7 rather than 10e-8 when it converges on the final step.

Overall Halley definitely converges faster. However there are problems.

The first is numerical instability. In the Newton code I used a value of dx=0.001. It worked great. In this code the second derivative requires dividing by dx squared, which would be 0.000001. Thats small. In fact it didn't work - I had to increase dx to 0.005 as a minimum otherwise the whole thing went crazy. Larger values of dx mean less accuracy to our derivatives, and things start to go wrong.

The second is that we're relying on more maths. Both Newton and Halley rely on the functions we're testing to be "well behaved", but there's just more to go wrong with that for Halley than Newton. There's also more code to go wrong - we're approximating more things, and that means more errors can creep in. In both cases there are a lot of things that need to go right to produce versions of the code that work reliably (the code here is for experimenting with and isn't designed to handle when things go wrong), but again theres more to go wrong in the more complex code.

The real problem however is that Halley converges faster PER STEP. But look at the work we're doing per step. It evaluates the function three times, rather than two, and then does a lot more maths on those values. Overall, although Halley uses less steps it probably uses about the same, or more cpu time. Computers are really good at doing lots of simple steps, so giving them the simpler code, and letting them run it fast, with lots of iterations is probably the better solution in the real world.

Even if this doesn't necessarily produce better results, it was still fun to try it, and there are a lot of good things to be learnt about implementing numerical analysis in code. We now return to our normal programming.

Solve any equation - Newton Rahpson

We're on a bit of a role of implementing advanced maths in Sniff, so lets keep going. So far we've looked at ways of calculating Sin and Square Roots, but perhaps we can write some code to solve any equation. Well maybe not any - sometimes we talk about "well behaved" functions, which means there some gotcha's that we're going to pretend don't happen (they're a bit like "smooth" surfaces in physics), but generally it will work for a lot of things we'd like to solve. This would also work just as well in Scratch, because Scratch is actually a lot more powerful than it looks.

Lets say we want to calculate a square root of 2. We can rewrite that a bit to get x*x=2, or x*x-2=0. We can generally rearrange any function to make it something on the left, equal to zero on the right. We're going to find f(x)=0, and to use another handy phrase, we can do this "without loss of generality" (this means it looks like we're simplifying things, but actually we're not - any equation can be rewritten in this form).

Here's some unknown, typical function. It crosses the X Axis, so there's a value x, where f(x)=0, and we'd like to find it. We've no idea where it crosses, so lets guess! What should we guess? Well for complex functions we maybe have to make a good guess, but for now lets just pick anything, and call our guess X0. We can find f(X0) and then (surprise) we realise that f(X0)!=0. d'oh!

So what next? Well remember how for a sin function we approximated it with a series of powers? Lets take that to its simplest form and just say that sin(x)=x. Guess what - that works if we're close to 0. In fact we can approximate any function with a straight line! After all straight lines are much easer.

If we replace f(x) with a straight line y=mx+c based around the f(X0), then we'd get a triangle like the one in the diagram. Now instead of solving our tricky function, we can solve our dumb approximation (cause its easier!): where does mx+c=0?

Well if we can find the gradient of the line: m, then we can fairly easily show that the the straight line will hit the X axis at:

X1=X0-(f(X0)/m).

To find the gradient to use for the line, so it's a good approximation to the curve at X0 we need to consider a tiny triangle: we know y at X0. Now move a little to the right (a distance we call delta x, or simply dx) and evaluate f(X0+dx). The gradient of the line we need is just (f(X0+dx)-f(X0))/dx. In other words how much we've changed in y divided be how far we've moved in x. The exact value of dx doesn't really matter - but ideally we'd like it to be as small as possible.

If you've done calculus you should have spotted that m is an approximation to f'(x), and we can often calculate f'(x) directly instead of this triangle/gradient stuff [strictly speaking using this approximation makes this the secant method rather than Newton's method, but that's probably not important]. However if you've not done calculus then just know that mathematicians sometimes write m as f'(x).

So X1 is the solution? Well sort of, except its the solution to the straight line, which is only an approximation to f(x), but just like in our sqrt code from the previous post, X1 is (probably) a better solution than X0. In that case, all we had to do we use each guess to calculate a better guess, and after a few goes round we get a pretty good answer. This gives us the final equation:

Given a guess Xn we can find a better guess by plugging it in to this equation. Intuitively we can see that if we get to the right answer then f(Xn)=0, so the value stops changing. This is called Newton's method, or Newton-Rahpson.

That was a lot of maths and no code! Lets write that in Sniff:

make x number

make f number

make dx number

make fx number

make m number

when start

.set x to 45

.set dx to 0.0001

.repeat 10

..broadcast calcF and wait

..set fx to f

..say join join [x] ":" [fx]

..

..change x by dx

..broadcast calcF and wait

..change x by -dx

..set m to (f-fx)/dx

..set x to x-fx/m

We start by setting our first guess x to 45, and loop around 10 times. We calculate the value of the function f, and store it in fx. We also print out x and fx so we can see how we're progressing.

To get a better guess, we increase x by dx and find f at this point (we then change x back, so we're back in the right place!). Now we can use the new value of f, along with the previous fx, to find the gradient m. Now we know everything we need to know to apply the newton-rahpson formula, get a better guess, and go back around the loop again.

So lets use that code to find the square root of 2. Remember the approach finds f(x)=0, so we write the problem as x*x-2=0. We can just drop that in to a script to calculate F, and let the rest of the code run:

when calcF

.set f to x*x-2

45:2023

22.2981:495.203

11.1071:121.369

5.65919:30.0264

3.00894:7.05371

1.83491:1.36691

1.46238:0.138568

1.41502:0.00227356

1.41421:7.15256e-07

1.41421:-1.19209e-07

Well the first thing we learn is that 45 is a pretty terrible guess for the square root of 2,  but after about 5 times around, replacing our current guess with a better one each time, things start to look promising - we've figured out its 1 and a bit! After nine times around we get the answer 1.41421 which is the correct answer to 5 decimal places. If you look at the 10th time around, and answer is the same when its printed out, but the error is even less.

But we didn't really think about what was going to be in calcF when we wrote the code, so we can drop anything we want in. What about cos:

when calcF

.set f to cos of x

45:0.707107

104.316:-0.247277

89.4999:0.0087276

90.004:-7.04406e-05

90:5.32632e-07

90:6.12323e-17

90:6.12323e-17

Well 45 is still a pretty terrible guess, but after only 4 times around the loop we've got the right answer and the error is effectively zero.

But we knew those, or could calculate them easily. Lets try something we don't know:

This is cubic, so its not so easy to solve. In fact equations like this are exactly the sort of thing Newton Rhapson is used for.  Lets just throw it at the code and see what happens:

45:95129

29.5867:27619.5

19.4132:8049.59

12.7117:2363.52

8.29118:698.161

5.37315:206.495

3.45087:60.4607

2.20288:17.1924

1.43379:4.62523

1.00966:1.05843

0.835976:0.145962

0.803114:0.00487411

0.801939:7.15256e-06

0.801938:2.38419e-07

0.801938:-1.19209e-07

Apparently it has a value of zero at about 0.8. Lets plot it and see (OS X has a neat graph plotting app in the Applications/Utilities folder):

Well so it does!

But wait a minute - that graph crosses the X Axis THREE times. There are three possible correct answers, but we've just found one. It all depends on that first guess (its also possible to make a really bad first guess, but we'll ignore that for now). It looks like there are crossing around -0.5 and -2.5, so lets try starting points of -1 and -2:

-1:1

-0.500083:-0.124813

-0.555568:0.00140119

-0.554958:-5.96046e-08

-0.554958:-5.96046e-08

-0.554958:-5.96046e-08

-2:1

-2.33341:-0.481945

-2.2531:-0.0317504

-2.247:-0.00013113

-2.24698:7.15256e-07

-2.24698:-9.53674e-07

-2.24698:7.15256e-07

-2.24698:-9.53674e-07

If we pick a value for out first guess that is close to a correct answer then it converges pretty quickly. If we plot a graph then we can usually use that to make a good first guess and from there the code does the rest!

Sniff (or Scratch) aren't my first choices of language for doing this kind of work - if you're serious about maths then you should be looking at something like Mathematica (which is FREE on raspberry pi, and justifies the code of a pi by itself), but its up to the job, and may be appropriate if you're teaching someone with strong maths, who hasn't learnt a more advanced programming language yet.

I think one of the main takeaways from this series of posts is how coding can really help in maths classes. Newton Rhapson and Taylor approximations are pretty abstract concepts, but they're ideally suited to being coded. Once implemented on code they become real things, rather than just equations on a page. 

Taylor Approximations - further maths!

Recently on the blog I've been looking at some neat mathematical algorithms that we can implement in Sniff. The code to calculate Pi and find square roots are probably beyond Sniff's target audience in terms of stuff you'd expect them to know, but the code is simple, and uses maths that's simple enough that it should be reasonable for a strong KS4 student to grasp whats going on if you choose to use them as exercises.

Today I'm going to do something that's a part of the A level further maths syllabus - Taylor approximations. These are definitely more advanced, as they're derived from calculus. To understand them you need to know about derivatives, and even to use them you need radians so they're not really going to be much use for younger kids. However they're something that older students should absolutely be encouraged to implement in code: they're tricky to calculate by hand, somewhat confusing, but yield great results when you actually run them on a computer.

I'm not going to explain the maths behind these, but will skip to the take home:

We can approximate complex functions like sin() with polynomials. In the case of sin() it works best when x is a small number (close to 0), but if we keep adding terms then the approximation just keeps getting better and better.

To implement that in Sniff there are a couple of gotcha's: there's no power function, and no factorial function! This makes the implementation more difficult , but actually forces us to think about the problem and produces a better final solution.

Lets start with a chat about factorials. Every programming class teaches you to implement factorials as:

def f(x)
    if x=1
        return 1
    return x*f(x-1)

clever huh? we define factorial of X based on factorial of X-1. Eventually X gets down to 1, which we already know the answer to (hint: its 1! by which I mean 1=1! by which I mean 1!=1! I mean!!!!! you get the idea!). This is recursion and its a big thing in programming classes.

Well not so fast... Scratch and Sniff can't do recursion so lets look at another way to write that:

make f number
when factorial
.make i number
.set f to 1
.repeat x using i
..set f to f * i

Here we've written the factorial without using recursion. Now its not "clever" - its just the obvious way of calculating 1x2x3x4x5x... and so on. If you've never been taught how to write a factorial, this is probably what you would probably write. Not only is it clearer, and simpler, its better! Each time the function f is called it creates a copy of x, and a bunch of other data. To find f(100) the recursive function calls itself  100 times, making 100 copies. Making those copies is slow, and uses up a lot of memory.

So lets get this clear: the fancy version that I was taught in programming classes is slower, uses more memory, and is harder to understand than the regular and "bleedin' obvious" version. Yep, sure is! The recursive version has NO advantages - the regular version is better in every respect. Advocates of recursion will point out that a really good compiler can maybe claw back the performance penalties, but that's really not the point.

The thing is recursion is a powerful idea. It's important, and does have uses (as an undergraduate I remember making fun of the physics students struggling to implement dozens of lines of code in Fortran that could be written elegantly in 6 lines of Pascal using recursion), but its also tricky to understand. Teachers in programming classes use factorial as an example because its simple, and they can't think of a good example of when you should use recursion, so they use a bad example. Kids: if you're ever in a programming class and get "taught" to write a factorial, ask why you can't use the simpler and faster version!

But lets get back to that sin function - in fact its going to be better not to break out a factorial function at all. Look again at that equation. For each term the the denomentator factorial increase by 2. The fastest way to calculate the denominator for a term is to perform two multiples on the previous denominator - much quicker and easier than than doing a full factorial every time.

Once we spot that we release that power of the numerator increases by two so we can get the next numerator by multiplying the current one by x*x - much cheaper than a power.

If we put that together we get:

make x number 
make s number 
make s number
when approxSin
.make numerator number
.make denominator number
.make count number
.set s to x
.set  numerator to x
.set denominator to 1
.set count to 1
.repeat 2
..set numerator to -numerator*x*x
..set denominator to denominator*(count+1)*(count+2)
..change count by 2
..
..change s by numerator/denominator

As this is an "advanced" example - at least mathematically, I've used local variables. These are optional, but once you realise you need them they're a big help - by defining variables inside a script we don't have to worry about other code using them at the same time - each script gets its own version.

The numerator starts of as x and the denominator as 1. Each time around the loop we multiply the numerator by -x*x. This not only increases the power by two, but flips the sign.

For the denominator we use the variable count to tell us how far the sequence of factorials we are. count starts as 1, then to get to the second term we multiply by 2 and 3. We increment count by 2 to indicate that the denominator is now 3!

Having calculated the new numerator and denominator we just add the fraction to the running total

And that's it. We go around a few times, and get the right answer. I've set repeat to 2 here (one less than the equation) as this is enough to get  useful answer, but leaves enough error to see that there's something interesting going on.

when start

.make count number

.repeat 10 using count

..set x to (pi/2)*0.1 *count

..broadcast approxSin and wait 

..say join "x"[x] 

..say join "approx:" [s]  

..say join "sin(x):"  [sin of (x/pi*180)]

..say join "error:" [s- sin of (x/pi*180)] 

..say ""

This script just tests it. We're working in radians, so set x to count from 0 to pi/2 in to 10 steps, calculating the approximation, the correct answer and the error. Note we have to use degrees with the built in Sniff sin function.

x0.15708

approx:0.156434

sin(x):0.156434

error:-1.3411e-07

x0.314159

approx:0.309017

sin(x):0.309017

error:-1.78814e-07

x0.628318

approx:0.587792

sin(x):0.587785

error:7.21216e-06

x0.942477

approx:0.809146

sin(x):0.809017

error:0.000128925

x1.25664

approx:0.952017

sin(x):0.951057

error:0.000960231

x1.5708

approx:1.00452

sin(x):1

error:0.00452483

Pretty impressive! The results are more accurate with small values of x (due to the way the approximation is derived), but even at pi/2 (90 degrees) the error is only 0.004. If we increase the number of iterations by only 1 additional term we can reduce that to 0.0001 which is probably good enough for a lot of applications (in fact Sniff uses code very similar to this to implement trig functions on the Parallax Propeller - the code uses much less memory than the official sin functions).

It's worth noting there's nothing here that couldn't be done just as easily in Scratch. I've used repeat using, and local variables which Scratch 1.4 doesn't have but they're purely for convenience, and the code would work just as well without them. So next time someone says that Scratch is too simplistic to do real work, as them if they've done Taylor approximations (or Fourier transforms, or Monte Carlo Simulation)