 Hey guys, welcome back skitstone series episode 27 topic today, you know, I'm always honest about this a filler topic complete waste of time Highly recommend you skip the video, but should be a quick one either way numerical integration or go through two different ways to implement this one is the trapezoidal method the other is Gaussian quadrature both pretty cool in their own right, but a bit of a speedrun for me in this video So let's get into it Integration, what is it basically? It's finding the area between a function and the x-axis or between functions perhaps And usually to do it so all of the speaking you have to do these dumb Math tricks and manipulations to make it work and as a result there has been no productive use case found for integration as Of the current month. However, if they find one which they won't I will amend the video in the description below So check that out if there's any changes But yeah, basically you get the signed area between this pink line and the white line In red here, and this is how you you draw that out, you know, and you go from a to b f of x dx Now that's great, but how can you enslave a computer to do that analysis for you? Couple options by the way, there's more than just this but here are my five ideas Option one. Why don't we render that function as a bitmap and then count the pixels or as a svg And count the pixels and then use that and then relate that to the area Good idea option to better idea How about we render as a bitmap and then print it out on my printer and then get my toner Capacity decrease and then relate that to the yeah, that's a good idea I would do that my printers broken though. Unfortunately option three. We could just guess Honestly, zero is a pretty good guess if you ask me, right? Area could be negative infinity or it could be positive infinity Zero is dead center, right? So not a not a bad guess if you ask me. So just return zero. That's fine option four If you really wanted to you could have probably approximate that area as You know as some set of geometric primitives that you can get the area of individually Maybe this is a triangle over here. Maybe this is a circle and that's a square This is some weird octagon looking thing, right? Some trapezoid going on here. Maybe another square Maybe a hot air balloon shape over here. I don't know Yeah, that would work if you get the areas of each one add them up here all set An option five Maybe there's some fancy math tricks that we could work with as well So we'll stick with option four and five if my printer worked would add in number two, but unfortunately here That's where we have so trapezoid method What's the idea you're going to estimate the area of the function Between the function of the x-axis using trapezoids in this case I have this really advanced function here y equals x squared and see we've got five trapezoids or really four trapezoids in One triangle at a one And so you get the area of each you add them up and you're done. So how does this work? Well for trapezoid it's basically a rectangle on on drugs So you get the midpoint or the average side length between left and right That's this one here Multiply that by the base like a rectangle and then you're done So this part of it gets the average side length. This is the base, you know easy peasy And so you do that you add them up and you're done in this case You do that for these five trapezoids yet four two point five the exact answer if you trust these mathematicians is 41.67 exactly I'm sure there's no more digits around here, but we're pretty close not that close but for five Actually six function calls not half bad to be you know a couple percent off like this And so we implement this in the trapezoid method Assembly code that we'll talk about in this video at the end But yeah, that's pretty straightforward. There are some opportunities to make this more efficient if you think about it, right? this side That's the right hand side of a four But the left hand side of a five and so maybe you could change this some maybe I don't know to Account for the fact that you're reusing that side Yeah, maybe it could be I'll leave that as an exercise to the reader So what are some pros and cons of this so? Symbolic or analytical Integration you need some either some fancy Symbolic toolbox like they have in Matlab or in Mathematica or in Maple Which costs money unless you pirate the software, which you probably can do Or you have to tie up an mathematician keep moving your basement feed-in crackers But again that costs money the rope isn't cheap, right? And so it's expensive to do this The trapezoid method is much cheaper because you have to just call the function you call this function six times And you're done no cracker is needed. So good idea. What about number two? The symbolic integration method in the actual set of symbols You can't integrate nothing you need me to give you this But here's a secret about society. There are no functions in society. No one's gonna email you a function and say Ingrade this that's not real did They in real life. There's no functions. It's just data And so you may have a sensor data you want to integrate which you can do with trapezoid method but you can't you don't need a Expression to integrate that also you can do black boxes. So you had a black box to integrate You could do so So if you had some finite element analysis software that takes ten minutes to run and gives you a single quantity output You can integrate that numerically speaking with this method without having any analytical expressions So that's up, you know, and also advantage the only problem with the Estimates is that it doesn't give you an exact answer without having a lot of trapezoids You can see here even with five trapezoids. We were still not that close to the exact answer Or as you get the exact answer pretty quickly with a symbolic analytical integral Can we do better answers obviously? yes, and For that one way we can get to better for certain functions is Gaussian quadrature and I'm not gonna get into all the details here. You can check out the Wikipedia page but basically you convert the integral into a sum of function evaluations that are weighted by a certain with a set of weights and so for endpoints and end weights for your Estimate you can get an exact integral for polynomials of order two and minus one or less So for a cubic with just two points you can exactly get the area of a cubic polynomial integral So how does it work? Well, you have this integral here. So a negative one to one. Those are your bounds of the integration You sum up the weights times the function evaluations. So you have n weights wi and n weights or endpoints x of i So yeah, that's how that works and then the weights and the points they come from some Legendre polynomial stuff too hard for me So we're just going to preload and also it's too hard to compute every time if you think about it If every time you wanted an integral you have to you know do all this fancy Luzon Legendre polynomial stuff it would take forever and so by preloading all the data in a table you're able to very quickly Compute these estimates of the integral And so yeah, the weights on the points themselves We're going to preload them as a table and you may notice that the bounds of integration are from negative one to one. However We can obviously translate any bounds to those bounds just by Adjusting our function, right? So if our bounds are negative 10 to 10 we can adjust we can scale our inputs and scale the function Accordingly to get the bounds to negative one to one to work in this way And so here's how that works basically if you think about it Let's say your bounds are negative 10 to 10 that means your range is 20 whereas the range of this function here this integral here is 2 and so this Expression here that converts your range. So B minus a to the Gaussian quadrature range of 2 Similarly you have that here. So you're scaling both the input that input point as well as the output By that domain change which is easy to do and then this over here that's the shift So the midpoint of one a negative one is zero But the midpoint of a B is just a plus B over two and so you have to shift your Point by that. So basically think about it. This just converts the Gaussian points into your points in your domain and then this converts the range Yeah, the the resultant value back to the size that you want it to be So, yeah, that's how that works and here's some more example more details So here are the points and weights for n equals one to five and so for a single point With evaluating the function at x equals zero with a weight of two that gives you the estimate of your Integral and that would be valid for two a minus one and being one. So that's one, right? This gives you the exact value of a linear integral and Then this would give you a cubic and so on and of course n equals five This this n would give you the ability to estimate with five points and five weights Which are all titulated here a polynomial of order nine Which is way more than you'd ever need because most things in real life far Cubic at the highest in my experience maybe quadratic most of the time cubic sometimes maybe Cortic one million years So, yeah, that's pretty good And so you can see here. Here's a random cubic polynomial 7x cubed whatever the Trapezoid method with two points right the two end points are being evaluated You can see that gives you an absolutely horrible estimate of the area, right? If this is the zero line, this is not even close to the area of the Polynomial integral, right? However, the Gaussian one with two points these two random points at negative root one over three positive root one over three You get those two values of the function Multiply those values each by one in this case add them together you get the exact integral It's pretty powerful stuff, right and so yeah, it's basically a cheat code that you can use to evaluate the integral of polynomials of You know decent order and of course if you wanted to you could extend this list or you could compute these values Formulatively right you can get the if you really wanted to you could get the points and weights for polynomial You know or n equals like a hundred, but I'm not sure why you would but you could do that You could go through and figure out the the weights and the the x i's For any order right if you really wanted to We're not gonna do that though So what is What are our examples today? We have two one for each method. Let's get into it So here you can see those two examples Let's go into the actual code first lib math Integration and so here is the trapezoidal method. So what do I have set up here? So basically? It's a function you pass in an RDI the function pointer. So we're going to pass in Some polynomial or some some other function that you want to evaluate and that function The only thing about that function is it has to take in an input in X and M zero and it returns this output Also in X and M zero so let's see you had y equals X squared as your function that you want to integrate or estimate the integral of You that function has to be set up in a way that it doesn't affect any other registers But X and M zero so you have to push them out to the stack pop them all off afterwards But yeah, it takes in X in X and M zero Squares it and returns that back in X and M zero. That's how that would work But then you also have to pass in in RSI the number of steps So again five steps was a decent estimate But a hundred would have been better, right? And so you pass in the number of steps you want to do number of trapezoids in RSI and then you pass in the bounce So X and M zero is the lower bound and one is the upper bound and that will estimate the area for you No problem How does this work? Well first things first. I push all the registers to the stack that I'm not going to affect I don't want to affect it and then I put them all back at the end down here and then what happens So basically you can see here. We are competing the step size So we're taking the upper and lower bound dividing that by the number of steps You see we convert from a integer to a double with this Instruction here so we divide by that giving our step size and X and M two then we have a running sum that we track X and M four and then we basically call left-hand side we evaluate left-hand side of our trapezoid and then We got with the right-hand side of the trapezoid and then we constantly recompute the area of each trapezoid And then we shift the right-hand side to the left-hand side for the next one and just add up the sums as we're going and so A very simple process. I'm sure you could figure out a similar or perhaps better way to do this So that's that Now what about the quadrature and so in this case? It's in many ways easier because there's less things to do, right? You don't have to do a hundred trapezoids only have to do two function calls right for order to estimate But it's a little bit different in how it works. So in this case same inputs pretty much pass in the Function pointer in RDI the same function has to you know have inputs and outputs in X M zero only Pass in the order in RSI So you want orders order for you pass in four for RSI and then the two bounds in X M zero X M M one and again It returns the area estimate in X M M zero and this works for polynomials But it can also give you a not half bad estimate for other things as well But I would use this exclusively polynomials if you could So again first things first you save the registers to the stack at the end you pop them off back But I will say while I'm down here. I'll show you what's going on. So we have pre-loaded You know before compile time literally in the binary all of the points and weights for orders one two three four and five So they're all saved here as 20 point numbers, right? You can see here's your plus and minus square root of three right here's your weights of one for the order order two estimate The challenging thing is now that although we're not Responsible for computing these values with the whatever we're doing to you know get the formulas to work We do have to you Actually figure out where these are in the table, right? So you you pass in an RSI the order that you want to use for the estimate order three Now you have to figure out where in my table is the order three points and weights How does that work? Well, not too bad But you do have to do some math to figure that out and so that's what this is doing right here It's basically figuring out how many bytes am I offset from? that Points and weights label this address in memory. How far are we from this? Like say I want to do order three polynomial order three estimate I should say how many bytes between this address and this address and this address and this address So can beat those two things and then at the end of that basically you have in rex the Pointer to the start of the weights and then you have in our DX the pointer to the start of the points So yeah, that's how that works then you just loop through evaluate the Polynomial or the function at those points multiply by the weight Compute the running sum here X of m4 and then and you're done right at the end We do have to shift so one thing you have to do though remember is you have to constantly convert back to the The range negative one to one in order to use this and so we can see here. There's some instructions that help you relate the Input X of i's into the function X of i's right convert the balance negative one to one to Your particular bounds X from zero to X from one and so there's some Function handling for that here But yeah, besides that it's pretty straightforward. You literally only have like this is your entire loop right from here to here That does all the the lifting for this algorithm So let's now look at the examples themselves. So Example a chap's way method. Let's just run it. So if I run this You can see that's the function I'm integrating at the top Y is negative 4x cube plus 3x squared plus 1 half x minus 51 and you can see Here are the estimates for the area with a different With a variable number of trapezoids for the estimate and so with one trapezoid We have a horrible estimate. It's the wrong sign completely of 240 and then as you get to higher amount trapezoids We get a better and better estimate and of course you don't have to Go through a loop like this you could start off at 25 you don't have to do one two three four in order However, one advantage to doing it in order like this you can get like kind of a convergence History or you can figure out. Well, hey, I was pretty here's my estimate for all these levels of trapezoids You know one two three four all the way to a hundred But you know we got close enough. I would say at 16 that was close enough for my liking And so from now on we're going to stop at 16 and so there is some benefits to going through and Doing a whole bunch of these with different step sizes, for example, but you don't have to do that So how does this work? Let's look at the code So Includes are very simple All we have is the trapezoid method as well as some printing stuff to see our output So very minimal and here's the function I was talking about before so this function has to take inputs And return outputs in x m zero but not touch any other registers And so here you can see this function just computes that polynomial expression And it puts the value back in x and zero Then how does this work? Well, we just loop through we print out at the beginning the integral itself Then we loop through print out some stuff run the method print out the answer And we do that from you can see r8 of one all the way to r8 of 25 And you can see for the trapezoid method function call We pass in the two bounds x1 zero x1 one put the function pointer and rdi and then the number of trapezoids in rsi So let's change this really quick to a bigger number. Let's do a thousand That's not thousand. That's a hundred And then run this you can see we've now approached negative 260 and that's the exact answer by the way But you can see we were pretty close at 25 But we were not there and actually we're still not there technically speaking, right? And we'll never exactly be there Unless by chance, but probably not for this particular polynomial So we're close, but we're not exactly there Now what about example b this one is the quadrature and remember this is a cubic function So we should be able to get this in I think about two minus one. So order two estimate would be capable of estimating an exact integral for order three polynomial So if I run this You can see the same polynomial for the You know function Order one wrong answer order two right answer three four and five all right answer The exact answer So if you think about this look what we have on the screen here for a thousand trapezoids a thousand Which is the thousand and one function calls best case, right? We're still not right with our estimate But with two function calls for the gaussian quadrature we get the exact answer And now that's true for up to a million you have a million trapezoids You'd still not be the exact answer. You'd be close be very close But with only two function calls you can get the exact answer. Is that not powerful? And this is not like You know, I didn't like put this in there that's literally Just But based off the numbers on the table the points and the weights get the exact Integral of 260 negative there. How's that working the kill? Let's take a look much the same thing It's pretty much just plug-and-play. We just changed the integration method from trapezoids to quadrature Oops, I had to fix that. That's bad documentation there Quadrature Again same polynomial The same printing only difference being is that you're not passing in number trapezoids. You're passing in the order Of the polynomial. We've only encoded the weights and the points up to order Five you could expand that you want to do up to 10 20 30 be my guest Um, but yeah, we only loop from r8 of 1 down to r8 of 5 but If you want to add more be my guest But yeah, so pretty cool if you ask me you're able to Get an exact integral of a polynomial expression with very few function calls Using gaussian quadrature So and by the way, a lot of things in the real world are either polynomials or can be very closely approximated as such And so there's a lot of applications for this if you're in the business of computing integrals, which is a waste of time Let's be honest. Um, but yeah, thanks for watching. Have a great day. See you in the next video