I've often thought that you could give people a huge leg up on calculus if you could first explain to them the discrete calculus, possibly with a programming language like Haskell. It would look like this:
We're going to start thinking of infinite sequences. A sequence starts with a value, and then there is another value, and so on, and so on. We might number these sequences and then write down rules for them; for example we might write $x_n = (n + 1) ^2$ to denote the sequence of square numbers starting with 1, `map (^2) [1..]`; if you're not familiar with these then we'll write out the first five as:
squares = 1 : 4 : 9 : 16 : 25 : map (^2) [6..]
I'm going to teach you one list-operator, it's a function which subtracts each number in a sequence from the one that comes before it in that sequence. For the first element, it just subtracts 0, so it doesn't do anything. This is written in Haskell as:
delta seq = zipWith (-) seq (0:seq)
When we use the above we get:
Prelude> take 30 (delta squares)
[1,3,5,7,9,11,13,15,17,19,21,23,25,27,29,31,33,35,37,39,41,43,45,47,49,51,53,55,57,59]
(Here the `take 30` just stops the computer from trying to give us all infinity of the results!)
Interestingly, the difference between square numbers is the set of odd numbers. This is not too hard to understand after you know it -- $(n + 1)^2 = n^2 + 2n + 1$, so of course after you subtract $n^2$ you get $2n + 1$.
Similarly when we start looking at the various counting sequences like [1..] (which again if you don't know Haskell is the same as 1:2:3:4:5:[6..], the infinite sequence of successive integers starting with 1) and [2..] we find:
Prelude> take 30 (delta [1..])
[1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1]
Prelude> take 30 (delta [2..])
[2,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1]
Prelude> take 30 (delta [3..])
[3,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1]
We notice that they're all the same except for that very first element, which basically just tells us the first value. But actually, `delta` has a special property which is that it is
invertible: just like you can invert `\x -> x - 3` with `\x -> x + 3`, you can invert this `zipWith (-)` beast with `scanl1 (+)`
Prelude> let sigma seq = scanl1 (+) seq
Prelude> take 30 (sigma (delta [3..]))
[3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32]
Prelude> take 20 (sigma (delta squares))
[1,4,9,16,25,36,49,64,81,100,121,144,169,196,225,256,289,324,361,400]
What is this magical beast `sigma`? What is `scanl1`? This is a function which processes a sequence by passing on an "accumulating value", but unlike `foldr` and `foldl` where you get that accumulating value at the end of the sequence, `scanl` gives it to you at each step of the journey, including the first. So like with foldl1, we start with the more general case of foldl or scanl:
scanl fn val list = val : case list of
[] -> []
x : xs -> scanl fn (fn val x) xs
So we always emit a sequence whose next element is `val` and whose follow-up sequence is either empty (if the source sequence is empty) or consists of scanl'ing the rest of the sequence with the same function but with val updated to `fn val x`. Each new element gets combined into the existing value with `fn` and the results are getting steadily emitted as it goes.
Then the "1" in `scanl1` means, just like the "1" in `foldl1`, "read the first value from the first element of the list":
scanl1 fn list = case list of
[] -> []
x : xs -> scanl fn x xs
So when we defined `sigma` as `scanl1 (+)` we just specified this `fn` to be one which adds the two values together. The above recomputation of the squares from the odd numbers actually just works like this:
1 = 1
1 + 3 = 1 + 5 = 4
1 + 3 + 5 = 4 + 5 = 9
1 + 3 + 5 + 7 = 9 + 7 = 16
1 + 3 + 5 + 7 + 9 = 16 + 9 = 25
1 + 3 + 5 + 7 + 9 + 11 = 25 + 11 = 36
Notice as the second column here indicates, we don't need to recompute all of these pluses each time: we can just hold onto the last number we've computed, the previous item from the last column, and add that to the new item to get the new value.
I'm asking you to notice the last point because if you study it carefully enough, that we generate each new value by adding the sequence to the last value, you'll start to ask, "wait, is `delta . sigma` also an identity transformation, just like we saw `sigma . delta` is?" And the answer is a resounding yes. This is actually a very common property for inverses, usually after you've proven that `f . g = id` then there's also an argument that `g . f = id`.
So, we've learned that just like `(+)` undoes `(-)` and vice versa, summing up a bunch of terms from some start point to a current place undoes taking differences between nearest-neighbors, and vice versa.
I haven't actually proven either side of this, but it's not too difficult: to see that delta undoes sigma, notice that for n > 0 we have $$
\sum_{k=0}^{n} x_k - \sum_{k=0}^{n-1} x_k = x_n
$$ (and the x_0 case is just x_0 - 0 = x_0, so that works too). To see that sigma undoes delta is a little more complicated, but just
write out the n'th term to see:$$
(x_0 - 0) + (x_1 - x_0) + (x_2 - x_1) + \dots + (x_{n-1} - x_{n-2}) + (x_n - x_{n-1}).
$$ Notice that except for the first and the last term, each term contains two parts: a +x_k which will cancel with the next term and a -x_{k-1} which will cancel with the previous term. This is called a "telescoping" series because it's really long but it collapses like a telescope due to perfect cancellations between adjacent terms. The only thing that's left at the end is $(-0) + x_n,$ which of course is just $x_n.$ So `sigma` has indeed undone `delta`.
To follow up, first I would probably define:
instance Num n => Num [n] where
(+) = zipWith (+); (*) = zipWith (*); (-) = zipWith (-)
abs = map abs; negate = map negate; signum = map signum
fromInteger = repeat . fromInteger
k !* seq = map (k *) seq
so that we have sequence algebra and we can multiply them by constants; then I would introduce people to basic discrete-calculus results like `delta (f + g) = delta f + delta g` and `delta (k !* f) = k !* delta f` and `delta (f * g) = f * delta g + g * delta f + delta f * delta g`.
Finally comes the big reveal: derivatives, to difference a continuous function, treat that function as if it were a straight line when you zoom in close enough, and to do this we have to divide `delta f` by some $d x$ to make it work; similarly integrals, to be the proper inverses to derivatives, contain a `(dx * )` operation. So if you start to define a grid of size dx over the x's that you care about, u < x < ∞, then each function corresponds to a sequence $f_k = f(x_k) = f(u + k * dx).$ Then `map (/ dx) . delta` generates a version of `delta` whose results converge, but not (in general) on the trivial answer 0; rather it converges on the "instantaneous slope" at x. Meanwhile `map (dx *) . sigma` perfectly undoes this, and calculates the area under a curve from u to x. From this, observe that the product rule gets easier due to these limits and introduce the formal definitions (limits, largest element of a partition going to zero for integrals; simple "h" parameter going to zero for derivatives).