cancel
Showing results for
Did you mean:  ## Introduction

My recent foray into solving complex numerical/engineering problems with DAX such as Kaplan-Meier Survival Curves and Linear Interpolation led me to wonder what other kinds of problems I might be able to solve using Power BI and DAX. This inquiry led me to remember one of my favorite engineering classes taught by Dr. Michael J. Rider at Ohio Northern University (ONU), Numerical Methods, and one of my favorite numerical method techniques, Runge-Kutta.

What follows is a story of utter frustration, some success, but ultimate failure to achieve what I was truly looking to accomplish. However, it is often the case that we learn as much or more from failure than we learn from success.  So, I felt it was important to share what I learned in trying to implement Runge-Kutta in Power BI and DAX.

## Background

In numerical analysis, Runge-Kutta is a family of iterative methods for approximating the solutions of ordinary differential equations. Now, if you have ever taken courses in differential equations, and I took four quarters worth of them at ONU, I'm certain that you can understand the desire to not solve them through traditional, pure mathematical means if at all possible. Such desires are seemingly near universal and have led to the development of things like Laplace Transforms, which essentially turn differential equation problems into Algebra problems. However, given the time-honored engineering principle of "Eh...close enough", numerical approximation techniques like Runge-Kutta were invented as well. Specifically, the Runge-Kutta methods were developed around 1900 by the German mathematicians C. Runge and M. W. Kutta.

Now, while there are an entire family of Runge-Kutta methods, the most widely used method is known as the fourth order Runge Kutta method (RK4). This method involves the calculation of a set of four numbers for each iteration of a "step size". These four calculations are then used to compute an approximation of the value of an equation at some point in time.

Consider the following problem:

y' = f(t,y)

y(t0) = α

We can define h to be the time "step size" and ti = t0 + ih. We can then make the following calculations:

w0 = α

k1 = h * f(ti, wi)

k2 = h * f(ti + h/2, wi + k1/2)

k3 = h * f(ti + h/2, wi + k2/2)

k4 = h * f(ti + h, wi + k3)

wi+1 = wi + 1/6 * (k1 + 2 * k2 + 2 * k3 + k4)

Such that wi y(ti).

Simple enough, or so it would seem...

Let's consider a specific example:

y' = y−t^2 + 1

y(0) = 0.5

Let's say that we are interested in the value of y for 0 ≤ t ≤ 2. If we choose a step size, h = 0.5 then it will take four steps to solve, t0 = 0, t1 = 0.5, t2 = 1, t3 = 1.5, t4 = 2.

Step 0

t0 = 0

w0 = 0.5.

Step 1

t1 = 0.5

k1 = h * f(t0,w0) = 0.5 * f(0,0.5) = 0.75

k2 = h * f(t0 + h/2,w0 + k1/2) = 0.5 * f(0.25,0.875) = 0.90625

k3 = h * f(t0 + h/2,w0 + k2/2) = 0.5 * f(0.25,0.953125) = 0.9453125

k4 = h * f(t0 + h,w0 + k3) = 0.5 * f(0.5,1.4453125) = 1.09765625

w1 = w0 + (k1 + 2 * k2 + 2 * k3 + k4)/6 = 1.425130208333333

Step 2

t2 = 1

k1 = h * f(t1,w1) = 0.5 * f(0.5,1.425130208333333) = 1.087565104166667

k2 = h * f(t1 + h/2,w1 + k1/2) = 0.5 * f(0.75,1.968912760416667) = 1.203206380208333

k3 = h * f(t1 + h/2,w1 + k2/2) = 0.5 * f(0.75,2.0267333984375) = 1.23211669921875

k4 = h * f(t1 + h,w1 + k3) = 0.5 * f(1,2.657246907552083) = 1.328623453776042

w2 = w1 + (k1 + 2 * k2 + 2 * k3 + k4)/6 = 2.639602661132812

Step 3

t3 = 1.5

k1 = h * f(t2,w2) = 0.5 * f(1,2.639602661132812) = 1.319801330566406

k2 = h * f(t2 + h/2,w2 + k1/2) = 0.5 * f(1.25,3.299503326416016) = 1.368501663208008

k3 = h * f(t2 + h/2,w2 + k2/2) = 0.5 * f(1.25,3.323853492736816) = 1.380676746368408

k4 = h * f(t2 + h,w2 + k3) = 0.5 * f(1.5,4.020279407501221) = 1.385139703750610

w3 = w2 + (k1 + 2 * k2 + 2 * k3 + k4)/6 = 4.006818970044454

Step 4

t4 = 2

k1 = h * f(t3,w3) = 0.5 * f(1.5,4.006818970044454) = 1.378409485022227

k2 = h * f(t3 + h/2,w3 + k1/2) = 0.5 * f(1.75,4.696023712555567) = 1.316761856277783

k3 = h * f(t3 + h/2,w3 + k2/2) = 0.5 * f(1.75,4.665199898183346) = 1.301349949091673

k4 = h * f(t3 + h,w3 + k3) = 0.5f * (2,5.308168919136127) = 1.154084459568063

w4 = w3 + (k1 + 2 * k2 + 2 * k3 + k4)/6 = 5.301605229265987

So, as one can see, the calculations at each step depend upon the solutions to those calculations in the previous step. Specifically, the calculation of k1 depends upon the calculation of w at the previous step and the calculation of the next w also depends upon the value of the previous calculation of w.

## The Solution...Sort of...

OK, so the original goal of this was to define a single Measure, w, that I could plot at each point in time t. The issue with this is that DAX hates circular references, I mean HAAAAAAAAATES them! With a passion. Try as I might, I simply could not get a single column or measure that could reference its previous value without having DAX just flatly refuse to cooperate and whine about circular references.

So, the solution I came up with ultimately had a single measure that I could plot, but required a number of additional measures as well, one for as many steps as required by the chosen step size.

Here are the calculations to solve the same problem as above.

First, define a measure h as: h = .5

Second, create a table using the following formula: RK4 = GENERATESERIES(0,2,[h])

Next, create the following measure: m_w0 = .5

Then this series of measures:

`m_w1 = VAR k1t = 0VAR k1w = [m_w0]VAR k1 = [h]*(k1w - k1t^2 + 1)VAR k2t = k1t + [h]/2VAR k2w = k1w + k1/2VAR k2 = [h]*(k2w - k2t^2 + 1)VAR k3t = k1t + [h]/2VAR k3w = k1w + k2/2VAR k3 = [h]*(k3w - k3t^2 + 1)VAR k4t = k1t + [h]VAR k4w = k1w + k3VAR k4 = [h]*(k4w - k4t^2 + 1)VAR w = k1w + (k1 + 2*k2 + 2*k3 + k4)/6RETURN wm_w2 =VAR k1t = [h]VAR k1w = [m_w1]VAR k1 = [h]*(k1w - k1t^2 + 1)VAR k2t = k1t + [h]/2VAR k2w = k1w + k1/2VAR k2 = [h]*(k2w - k2t^2 + 1)VAR k3t = k1t + [h]/2VAR k3w = k1w + k2/2VAR k3 = [h]*(k3w - k3t^2 + 1)VAR k4t = k1t + [h]VAR k4w = k1w + k3VAR k4 = [h]*(k4w - k4t^2 + 1)VAR w = k1w + (k1 + 2*k2 + 2*k3 + k4)/6RETURN wm_w3 =VAR k1t = 2*[h]VAR k1w = [m_w2]VAR k1 = [h]*(k1w - k1t^2 + 1)VAR k2t = k1t + [h]/2VAR k2w = k1w + k1/2VAR k2 = [h]*(k2w - k2t^2 + 1)VAR k3t = k1t + [h]/2VAR k3w = k1w + k2/2VAR k3 = [h]*(k3w - k3t^2 + 1)VAR k4t = k1t + [h]VAR k4w = k1w + k3VAR k4 = [h]*(k4w - k4t^2 + 1)VAR w = k1w + (k1 + 2*k2 + 2*k3 + k4)/6RETURN wm_w4 =VAR k1t = 3*[h]VAR k1w = [m_w3]VAR k1 = [h]*(k1w - k1t^2 + 1)VAR k2t = k1t + [h]/2VAR k2w = k1w + k1/2VAR k2 = [h]*(k2w - k2t^2 + 1)VAR k3t = k1t + [h]/2VAR k3w = k1w + k2/2VAR k3 = [h]*(k3w - k3t^2 + 1)VAR k4t = k1t + [h]VAR k4w = k1w + k3VAR k4 = [h]*(k4w - k4t^2 + 1)VAR w = k1w + (k1 + 2*k2 + 2*k3 + k4)/6RETURN w`

Finally, create this last measure:

`m_w = SWITCH(MAX(RK4[t]),0,[m_w0],1*[h],[m_w1],2*[h],[m_w2],3*[h],[m_w3],4*[h],[m_w4])`

We can now plot our approximate values for our function by using RK4[t] and m_w: Conclusion

Now, this solution passes over the hours spent in hair pulling frustration trying to thwart DAX's circular reference checking. Trust me, DAX has the circular reference checking thing down to a science and is apparently hack-proof. So, while I was ultimately able to achieve the visualization, I do not consider this exercise a true success. The reason is that the solution is overly complicated and limited given what was trying to be accomplished. The fact that you have to create a separate measure for each step in the Runge-Kutta method severely limits the usefulness of the approach. What I REALLY wanted to be able to do would have been a single measure that used a mythical PREVIOUSVALUE DAX function like:

```m_w =
VAR k1t = [h]
VAR k1w = PREVIOUSVALUE(m_w)
VAR k1 = [h]*(k1w - k1t^2 + 1)
VAR k2t = k1t + [h]/2
VAR k2w = k1w + k1/2
VAR k2 = [h]*(k2w - k2t^2 + 1)
VAR k3t = k1t + [h]/
VAR k3 = [h]*(k3w - k3t^2 + 1)
VAR k4t = k1t + [h]
VAR k4w = k1w + k3
VAR k4 = [h]*(k4w - k4t^2 + 1)
VAR w = k1w + (k1 + 2*k2 + 2*k3 + k4)/6
RETURN w```

This would have allowed me to choose any step size and range of values and still have a single measure. Unfortunately, to the best of my knowledge, no such DAX function exists and therefore, perhaps this particular kind of problem should be left to the MatLab's of the world...