Arduino breathalizer (with IR remote!)

Ages ago I picked up one of these MQ3 alcohol sensors. They're part of a whole series of similar sensors which can detect the levels of various gases. Most of them are a few pounds so if you have an application where you might be worried about specific chemicals then you can pick an appropriate one. The exception to this is the Carbon Dioxide sensor which is hard to get and expensive, which is unfortunate as it would actually be interesting so see how levels change in rooms as they fill up with people. As I just wanted to play around with it, I picked the alcohol vapour sensor - I figured alcohol  was something I had easy access to that I could use safely.

Though you can get them separately, its easiest to buy the sensor already mounted on a board, like this one.  These generally have four pins. Two are for gnd/power, and the other two are outputs. Of those one is digital, in that it either turns on or off when the gas level reaches a certain threshold, setable by a small pot on the back. If you just want an alarm, you can hook this straight to a piezo buzzer, and you're good to go.

We're more interested in the other output - the analog output which produces a voltage dependant on the gas level. This will connect straight in to an Arduino, and as I just got a multifunction shield, I though I'd use that.

We'll keep all of the display code from the last two posts, and just write the new code:

make alcohol analog input A5

make min number

when start

.set buzzer to high #or it will make noise!

.forever

..set message to [alcohol*100]

..wait 0.1 secs

Easy! we just read the value (which will have a value between 0 and 1) and scale by 100.

But what does that value mean?

Others have conducted this experiment and failed to accurately calibrate these sensors. This is in part due to them forgetting to write down the results during the more intensive parts of the test regime, but there's general consensus that its hard to get any specific numbers.

In my tests the reading generally started around 20 when I turned it on, but as the sensor warmed up (its got a heating element in it), the value dropped to around 3. Putting it near the top of a bottle of vodka instantly spiked the value to around 20 again, but it dropped back again quickly when moved. Drinking a tiny amount (in the name of SCIENCE!) so there would be alcohol on my breath, produced similar results (15-20), but the level on my breath dropped back down again fairly quicky, as the vapour dissipated from my breath.

A more extended test involving half a bottle of red wine (no one said this was going to be easy) produced a reading of around 15 some time after I'd stopped drinking. As there was alcohol in my bloodstream, the level on my breath was fairly constant.

Putting these together I came up with a fairly simple ad-hoc approach which should detect the basics - record the minimum reading the sensor has ever returned as a baseline. If you blow double that you've been drinking but are probably quite sober. If you blow 4 or 5 times the minimum you're probably having a good time. If you blow 8 times the minimum, then you're wasted!

make min number

when start

.set buzzer to high #or it will make noise!

.set min to 1

.forever

..set message to [alcohol*100]

..if alcohol < min

...set min to alcohol

..set buzzer to not alcohol > 4*min

..wait 0.1 secs

That's pretty easy to code. Each time around we just check against the minimum, and update the minimum if the new version is lower. We check if the reading is higher than 4 times the minimum and sound the buzzer if it is. You can adjust the threshold for this depending on if you want to scare your kids, or set a frat house drinking challenge, but 4 seemed a fair value to me to call someone drunk and set of the alarm.

It's important to remember that this calibration is pretty arbitrary, and even if it was correct it only applies to the sensor I bought, not the one you have - under no circumstances you should you think you're safe to drive (or do pretty much anything) because some circuit you found online says so!

So given this is just for fun, lets have some fun with it!

make min number

make fakeResults boolean

when start

.set fakeResults to no

.set buzzer to high #or it will make noise!

.set min to 1

.forever

..if alcohol < min

...set min to alcohol

..

..if fakeResults

...set message to [alcohol*5*100]

...set buzzer to off

..else

...set message to [alcohol*100]

...set buzzer to not alcohol > 4*min

..wait 0.1 secs

Here I've added a fakeResults flag. When fakeResults is true it sets the reading to 5 times what it should be, and sounds the alarm (remember the buzzer is logic low). In regular mode it works as before.

So how do we control the mode? Well we could just use a button, but it might be a giveaway that you're cheating. The multifunction shield has a socket for an IR receiver. If you get a tsop4838 (about $1 on eBay) then it will plug straight in to the left three empty sockets. To use it all we need to do is make a receiveIR device on pin 2. Then tell it to "read". It sets a variable keyPressed to a different number for each key (just write a program to print out keyPressed to find the appropriate values for any remote).

make irSensor receiveIR device 2

make keyPressed number

when start

.forever

..tell irSensor to "read"

..if keyPressed=18615

...set fakeResults to not fakeResults

..wait 0.1 secs

Then we check if a key is pressed and use it to toggle fakeMode.

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.