Solving RK4 Integration with Matlab for Paths

[skyflashings]

Hi, I am trying to apply the fourth order Runge-Kutta integration method to a problem but I'm not sure exactly how to do it.

I have a matrix which represents a coordinate system (x, y) with varying degrees of height at each point. Using Matlab, I made a matrix with a gradient (u and v components) across the system.

Now, I want to place a particle at some point in the system and watch it's path be 'pushed' around by the gradients over time. You can imagine this like the current in an ocean pushing a boat around I guess.

I have applied Euler's method to this, and it was successful. Of course, the problem is Eulers is not that accurate, so I want to try RK4 to get a better path approximation. This is where I get lost. For Euler's, I used pathX(i) = pathX(i-1) + u(x, y) * dt and pathY(i) = pathY(i-1) + v(x, y) * dt. Since u and v are the partials (gradients) of the original height data, it seems to work fine.

For RK4, I need more than just the gradient I think. But I don't really know how to approach it. I am finding it difficult because I don't really know of a 'function' per se to plug values into. What does f(xi + (1/2)*h, yi + (1/2)*k1*h) mean in my context? Is it referring to solving the function using a different kind of stepped gradient or something else?

I would appreciate any help. If I am confusing or didn't explain the situation well, please ask! Thanks.

 

[I like Serena]

Hi skyflashings! Welcome to PF!

Your problem seems to be that you made a matrix of gradients.
That won't work with Runge Kutta.
You need functions u(x,y) and v(x,y) to find the slopes at an arbitrary point (x,y)

With Euler you only need 1 slope at a predictable point.

With Runge-Kutta you need 4 slopes at variable points.

What does f(xi + (1/2)*h, yi + (1/2)*k1*h) mean in my context?

xi is the time at iteration i
h is your time step dt
yi is the position vector (x,y) of your particle at iteration i
k1 is one of the 4 intermediary points where the gradient needs to be calculated (for iteration i)

So the function result of your f here is a gradient given by:
f.x = u(yi.x + (1/2)*k1.x*dt, yi.y + (1/2)*k1.y*dt)
f.y = v(yi.x + (1/2)*k1.x*dt, yi.y + (1/2)*k1.y*dt)

 

[skyflashings]

You need functions u(x,y) and v(x,y) to find the slopes at an arbitrary point (x,y)

Thank you so much for replying. I'm starting to feel like I'm getting a better grasp of it. However, I have a question I hope you can help me with.

For the functions u and v as you described, how can I find those slopes? All I have is a 2-d matrix with number values at each x,y position. I can kind of understand how a gradient is made since it is vectors pointing in the direction of the increasing values of the matrix. But how would I find the slope of that matrix at some point? I think I don't really understand how the gradient is done. Thanks again!

 

[I like Serena]

You wrote that you made a matrix with a gradient (u and v components) across the system.
How did you calculate it?

Since you apparently only have a matrix of heights, I guess the logical way would be to interpolate the height matrix with a cubic spline surface.

From the interpolated surface you could calculate a gradient at any (x,y) and use that for RK.

But I'm guessing here because I do not know what the context of your problem is.

 

[skyflashings]

Since you apparently only have a matrix of heights, I guess the logical way would be to interpolate the height matrix with a cubic spline surface.

From the interpolated surface you could calculate a gradient a any (x,y) and use that for RK.

I believe you have just pointed me in the right direction! Actually, the Matlab program created the gradient automatically for me. However, when looking around for interpolation methods in MATLAB I came across interp2, which will interpolate for you when you give it an original discrete matrix. So, I think that will indeed serve my purpose.

Just one more question if you don't mind: Should I take the gradient I already have and interpolate the in-between value I need? Or should I just interpolate on the original matrix and then perform a gradient on that? Sorry to ask so many questions! I really appreciate your help; I'm sure it would have taken me a much longer time to get this far.

 

[I like Serena]

Just one more question if you don't mind: Should I take the gradient I already have and interpolate the in-between value I need? Or should I just interpolate on the original matrix and then perform a gradient on that? Sorry to ask so many questions! I really appreciate your help; I'm sure it would have taken me a much longer time to get this far.

I understand you got your gradients from the heights matrix.
In that case you do not want to use it.
The gradients should be calculated from the bicubic spline interpolations, that in turn should be based on your input heights matrix.

Btw, I do not know MATLAB, so I can't really advise you there.
Preferably you would have a function that accepts an x and y and a matrix of heights, and yields a gradient (dx, dy) through bicubic spline interpolation.
I suspect however, that this is too advanced for MATLAB.

 

[skyflashings]

The gradients should be calculated from the bicubic spline interpolations, that int turn should be based on your input heights matrix.

I see. So perform the gradient on the interpolation, not the other way around. That certainly makes sense.

I really appreciate the direction, and your time in explaining it to me. I'm thoroughly impressed :)

 

[skyflashings]

If I interpolate first, I get a single height value as a result right? How can I apply a gradient across that value in relation to the rest of the matrix? Do I somehow add it into the original matrix and do a gradient from there?

 

[I like Serena]

If I interpolate first, I get a single height value as a result right? How can I apply a gradient across that value in relation to the rest of the matrix? Do I somehow add it into the original matrix and do a gradient from there?

I guess you can't do that.
You're limited by the tools MATLAB offers.

As an alternative you might let MATLAB make a matrix of gradients as before, and use that to spline-interpolate a gradient at a, by RK requested, point (x,y).

 

[skyflashings]

As an alternative you might let MATLAB make a matrix of gradients as before, and use that to spline-interpolate a gradient at a, by RK requested, point (x,y).

Ok, I'll try that approach, thanks again.

 

[skyflashings]

I just finished the implementation of taking the spline on the gradient matrix for RK4 and it works wonderfully. A very smooth path is made, but the path is still relatively similar to my Euler implementation. Unfortunately, I was hoping to get drastically different results.

Thank you so much for your time 'I Like Serena'. You've saved me quite a headache :)

 

[I like Serena]

I'm glad to hear that you got it to work!

Of course the accuracy is only as good as your gradient matrix and the interpolation thereof.

At least you should see that even with a bigger step size, RK4 has a way better accuracy than Euler, and its effective results are calculated much faster, despite the extra calculation overhead.

 

[skyflashings]

I'm glad to hear that you got it to work!

Of course the accuracy is only as good as your gradient matrix and the interpolation thereof.

At least you should see that even with a bigger step size, RK4 has a way better accuracy than Euler, and its effective results are calculated much faster, despite the extra calculation overhead.

Yes, I was quite surprised to see how quickly the result was computed!

In the path simulation, the particle spins around a loop in the current about three times and then exits off to the right. But if I use a smaller time step dt the particle stays to spin nearly 6 times or more before breaking out to the right. The end result is the same, but I find it interesting that the time step can influence RK4 so much. Is that usually the case?

 

[I like Serena]

Yes, I was quite surprised to see how quickly the result was computed!

In the path simulation, the particle spins around a loop in the current about three times and then exits off to the right. But if I use a smaller time step dt the particle stays to spin nearly 6 times or more before breaking out to the right. The end result is the same, but I find it interesting that the time step can influence RK4 so much. Is that usually the case?

Actually you should see that Euler is really bad at keeping a particle in a orbit that's supposed to be stable.
RK does a much better job and that with a much greater step size.
Of course, you'll still have an approximation error (4th order instead of 1st order).

The trick is too reduce your step size until the path it takes does not change noticeably any more.
And if you have sharp angles in your path, you might consider a variable step size that adjusts to the curvature of the trajectory.