# Fourth order Runge–Kutta in C# - two differential equations

1. Jul 20, 2011

### Itosu

Hello!

I am trying to solve two differential, non-linear equations in C# using fourth-order Runge-Kutha method.

I implemented them in C#:
Code (Text):

private const int ALFA = 7600;
private const int MU = 7200;
private int BETA = 1430;//int.Parse(HttpContext.Current.Session["BETA"].ToString());
private const int LAMBDA = 2470;
private int V = 139000;//int.Parse(HttpContext.Current.Session["V"].ToString());
private const double FI = 0.51;
private const double TETA = 2.5;
private const int QL = 8400;
private const double H = 0.01;

public double x = 0.81;
public double y = 0.055;
// ;
public ArrayList yList;
public ArrayList xList;

public void Start()
{
xList = new ArrayList();
yList = new ArrayList();
double t = 0;

//   xList.
// y = Y0;
for (double krok = 0; krok <= 500; krok++)
{
t = krok / 100;
x = rungeDX(x, t);
y = rungeDY(y, t);
//x = X0 + (i * H);
//y = runge(x, y);

}

}

public double doplyw(double t)
{
if (t >= 0.5 && t < 0.75)
return 100000;
else return 0;
}

private double dx(double x, double y, double t)
{
double result;
if (x > TETA)
result = (doplyw(t) + QL - LAMBDA * x - V * x * y - MU * (x - TETA)) / 15000;
else
result = (doplyw(t) + QL - LAMBDA * x - V * x * y) / 15000;
return result;
}

private double dy(double x, double y)
{
double result;
if (x > FI)
result = (-ALFA * y + BETA * (x - FI)) / 15000;
else
result = (-ALFA * y) / 15000;
return result;
}

double rungeDX(double x, double t)
{
double K1 = dx(this.x, y, t);
double K2 = dx((this.x + 0.5 * H), (y + 0.5 * K1), t);
double K3 = dx((this.x + 0.5 * H), (y + 0.5 * K2), t);
double K4 = dx((this.x + H), (y + K3), t);
double rest = (K1 + 2 * K2 + 2 * K3 + K4) * H;
double runge = this.x + rest;

return runge;
}

double rungeDY(double x, double t)
{
double K1 = dy(x, y);
double K2 = dy((x + 1 / 2 * H), (y + 1 / 2 * K1));
double K3 = dy((x + 1 / 2 * H), (y + 1 / 2 * K2));
double K4 = dy((x + H), (y + K3));
double rest = (K1 + 2 * K2 + 2 * K3 + K4) * H;
double runge = y + rest;
return runge;
}

}

However the results are not that that should be (in comparison to Simulink model).
Do you have any ideas?

#### Attached Files:

• ###### Capture5.png
File size:
3.8 KB
Views:
538
2. Jul 20, 2011

### I like Serena

Welcome to PF, Itosu!

You have mixed up to the function of x and t in your formulas.

Btw, it's a bit weird to have an argument x to a function if you're not going to use it, but use this.x instead.

Code (Text):

double rungeDX(double x, double t)
{
double K1 = dx(this.x, y, t);
double K2 = dx((this.x + 0.5 * H), (y + 0.5 * K1), t);

should be something like:
Code (Text):

double rungeDX(double x, double t)
{
double K1x = dx(this.x, y, t);
double K1y = dy(this.x, y, t);
double K2x = dx((this.x  + 0.5 * K1x), (y + 0.5 * K1y), t + 0.5 * H);

3. Jul 21, 2011

### Itosu

Last edited: Jul 21, 2011
4. Jul 21, 2011

### Itosu

The newest version is like this:
Code (Text):

public void rungeDX(double x, double y, double t)
{
double K1x = dx(x, y, t);
double K1y = dy(x, y);

double K2x = dx((x + 0.5 * H), (y + 0.5 * K1x * H), t + 0.5 * H);
double K2y = dy((x + 0.5 * H), (y + 0.5 * K1y*H));

double K3x = dx((x + 0.5 * H), (y + 0.5 * K2x * H), t + 0.5 * H);
double K3y = dy((x + 0.5 * H), (y + 0.5 * K2y * H));

double K4x = dx((x + H), (y + K3x * H), t + 0.5 * H);
double K4y = dy((x + H), (y + K3y * H));

double restX = (K1x + 2 * K2x + 2 * K3x + K4x) * H / 6;
double restY = (K1y + 2 * K2y + 2 * K3y + K4y) * H / 6;

this.x = x + restX;
this.y = y + restY;

}

and looks like the one in Simulink :)

Last edited: Jul 21, 2011
5. Jul 21, 2011

### I like Serena

Yes, I have some comments. :)
It should be like this:

Code (Text):

public void rungeDX(double x, double y, double t)
{
double K1x = dx(x, y, t);
double K1y = dy(x, y);

double K2x = dx((x + 0.5 * K1x), (y + 0.5 * K1y), t + 0.5 * H);
double K2y = dy((x + 0.5 * K1x), (y + 0.5 * K1y));

double K3x = dx((x + 0.5 * K2x), (y + 0.5 * K2y), t + 0.5 * H);
double K3y = dy((x + 0.5 * K2x), (y + 0.5 * K2y));

double K4x = dx((x + K3x), (y + K3y), t + H);
double K4y = dy((x + K3x), (y + K3y));

double restX = (K1x + 2 * K2x + 2 * K3x + K4x) / 6;
double restY = (K1y + 2 * K2y + 2 * K3y + K4y) / 6;

this.x = x + restX;
this.y = y + restY;

}

Cheers!

6. Jul 21, 2011

### Itosu

Thank you for your response, but unfortunately I do not understand why :)
My solution gives me exactly the same graphs that a simulation in Simulink, and your's does not.
However, you are the master, not me, so I am afraid that you are correct, but please explain :)

7. Jul 21, 2011

### I like Serena

I checked the http://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods" [Broken]and I see we need another modification.
It should be:

Code (Text):

public void rungeDX(double x, double y, double t)
{
double K1x = H * dx(x, y, t);
double K1y = H * dy(x, y);

double K2x = H * dx((x + 0.5 * K1x), (y + 0.5 * K1y), t + 0.5 * H);
double K2y = H * dy((x + 0.5 * K1x), (y + 0.5 * K1y));

double K3x = H * dx((x + 0.5 * K2x), (y + 0.5 * K2y), t + 0.5 * H);
double K3y = H * dy((x + 0.5 * K2x), (y + 0.5 * K2y));

double K4x = H * dx((x + K3x), (y + K3y), t + H);
double K4y = H * dy((x + K3x), (y + K3y));

double restX = (K1x + 2 * K2x + 2 * K3x + K4x) / 6;
double restY = (K1y + 2 * K2y + 2 * K3y + K4y) / 6;

this.x = x + restX;
this.y = y + restY;
}

You can see the regular form on wikipedia, where they define the problem as:
y' = f(t, y)
and they use a time step h.

What you need to be aware of, is that this y is actually the vector (x,y).
So everywhere where you read "y" in the wikipedia article, you need to replace it by "x, y".
The intermediary constants k1, k2, k3, k4 are also vectors, meaning they have an x and an y component.
And instead of a function f(t, y) you have 2 functions dx(x, y, t) and dy(x, y), which together form a vector again.

Furthermore in your program you have reversed the order of the position vector and the time.

Last edited by a moderator: May 5, 2017
8. Jul 21, 2011

### Itosu

You are totally correct. Thank you!
Do you know of any possibilities to calculate an error of such solution?

9. Jul 21, 2011

### I like Serena

You're welcome.

The typical method to estimate the error is the following.

1. Calculate the step with stepsize H yielding a new point (x1, y1)
2. Recalculate the step with stepsize H/2, and then calculate another step with H/2, yielding a point (x2, y2).
3. The difference between these 2 points is the estimation of the error.

The resulting error is accurate up to O(h5).

This method is typically used for the adaptive step size Runge Kutta method.
In this method the step size is dynamically adjusted (halving it or doubling it) to keep the error within a preset bound.

10. Jul 22, 2011

### Itosu

Hi :)
Thank you again!

And one more question...
I would like to calculate Area Under the Curve that I printed after having resolved the equations. I know I can do this with rectangle method, adding the areas of each rectangle, however I have a premonition that it has huge error then. What do you think?
Is there another possibility to calculate AUC?

11. Jul 23, 2011

### I like Serena

Rectangle method would work if you take your step size small enough.
To increase the accuracy you need to weigh your function evaluations differently.

On wikipedia you can find a list of methods:
http://en.wikipedia.org/wiki/Newton%E2%80%93Cotes_formulas" [Broken]

You would want to use the composite rule, meaning you subdivide your interval in many small intervals, and combine that with one of the schemes listed.

If f0, f1, f2, ..., fn are your function values with a step size H between them.
You have for instance the options:

Trapezoid rule (error is O(H3) which is less accurate than Runge-Kutta).
n * H * (f0/2 + f1 + f2 + f3 + ... + f(n-1) + fn/2)

Simpson's rule with n even (error is O(H5) which is more accurate than Runge-Kutta):
(n * H / 6) * (f0 + 4 f1 + 2 f2 + 4 f3 + 2 f4 + ... + 2 f(n-1) + fn)

Of course you can make it as complex as you'd like, but I guess first you would need to find out what problems you run into exactly if any.
If you have discontinuities or wild fluctuations in your derivatives you may need to take special care.
As it is I, believe you have a discontinuity, so you may consider doing the integration in such a way that you first integrate up to the discontinuity if that's feasible, and then add a second integration starting from the discontinuity.

If necessary, there are adaptive integration methods that may help, but if you're programming this yourself, I would only go there if you have need.
As an alternative, you may consider using math packages.
There are bound to be a few available for C#.NET.

Last edited by a moderator: May 5, 2017
12. Jul 23, 2011

### AlephZero

Erm.... "calculating Area Under the Curve" is the same as "integrating".

You just got your code for integration working, so why do you want to use a different method for the next part of the problem?

13. Jul 24, 2011

### Itosu

What do you mean? Can you explain?

14. Jul 24, 2011

### I like Serena

You did not really explain what you need the area under "the curve" for.
Depending on what you want exactly, it's possible to let Runge-Kutta do the work for you.

What you have is functions dx and dy that define a vector field.
You're using RK to find a curve with a gradient that matches these functions.
Effectively you are integrating the functions dx and dy.

If I understand you correctly, you want to know the area under the curve you just found.
This means that you are integrating the function again.
You could set this up as a new second order differential equation, using RK to find the result, instead of the first order differential equation that you are solving now.

It depends on what your results are for, to decide what's best and easiest in your case.
If you really only want the curve, and only want to display the area under the curve as an informative extra, I would integrate as I described in my previous post.
If you actually want the area under the curve, you can set up a different calculation and let RK find that for you.

15. Jul 24, 2011

### Itosu

My results are for calculations - I need to know the exact value od AUC (not for printing, just for forward calculations).
Do you mean that if I use the same (4th order) RK solver I will get AUC? Could you provide more details?

Thank you very much!

16. Jul 24, 2011

### I like Serena

Okay, I've just worked out the math, so here goes.

If I've understood everything correctly, what you would need to do is introduce two extra variables in your class:
Y (to keep track of the surface) and
t (to keep track of time).
Furthermore you would have to rewrite the rungeDx step function (shall we call it the rungeStep function?), and adjust the Start function.
And that's all!

You may not understand everything I write below, but could you at least answer the questions set in bold?

@AlephZero: Feel free to comment.

First let's rename your functions dx and dy to fdx and fdy to avoid ambiguity.
So you would have functions fdx(x, y, t) and fdy(x, y) that calculate the time derivatives of respectively x and y at some point (x, y) at time t.

Your current problem is then defined by de differential equation (DE):
dx/dt = fdx(x, y, t)
dy/dt = fdy(x, y)​
Btw, the only dependency on t, is apparently an impulse function for 0.5 <= t < 0.75 that is added to fdx(x, y, 0).

Until now you used RK to find a solution (x(t), y(t)) given a starting point (x0, y0, t=0).

If I understand correctly, what you're actually looking for is the function Y(x) which I define to be the area under the curve y(x) between x=x0 and some value of x.
Can you confirm this?

To find Y(x) with RK we need to set up a new DE which is:
dY/dx = gdY(x, Y)​
So the problem becomes finding this function gdY, and also take care of the impulse dependency on t.

We can do this by hauling y back in again, which is the derivative of Y(x), and also getting t back:
dY/dx = y
dy/dx = gdy(x, y, Y, t)
dt/dx = gdt(x, y, Y, t)​

Differentiation of y(x) gives:
dy/dx = (dy/dt) / (dx/dt) = fdy(x, y) / fdx(x, y, t)

Your problem would be solved by feeding the following system to RK:
dY/dx = y
dy/dx = fdy(x, y) / fdx(x, y, t)
dt/dx = 1 / fdx(x, y, t)​

There are a few boundary conditions that have to be satisfied though.
1. t must be a well defined function of x.
2. y must be a well defined function of x
3. dx/dt must never be zero.

So here I'm getting back to you.
Can I assume that x strictly increases with time?
Because that would satisfy all these criteria.

For your program this would mean introducing an extra variable Y to keep track of the area, and an extra variable t to keep track of time.
Rewriting the Runge step function, so it steps through x instead of through t.
Furthermore calculate successive values of (y, Y, t) with increasing x, instead of (x, y) with increasing t.

I can do that for you, but before I do, I'd like to know if this is useful for you, and if this is what you want.
Is it?

17. Jul 24, 2011

### Itosu

Yes.

Sample graph is like the one in the attachment. So, as you see x increases, however after some time it dramatically increases, because what I am (We are) creating is glucose-insuline model, and when glucose level increases more insuline is produced and it makes glucose level decrease.

AUC is very useful for me because it helps to model the effect of glycemic index, so if you had a little time I would be very grateful for the calculations. I need more time to think over what you just wrote :)

#### Attached Files:

• ###### AUCref.PNG
File size:
5.8 KB
Views:
417
18. Jul 24, 2011

### I like Serena

What do the axes in your graph represent?
Is it an x-t-diagram or what?

And what do your variables x, y, and Y represent?
I presume they are related to glucose level and/or insuline level?

I presume "Czas" is time in hours?
Is "Glikemia" glucose level?

19. Jul 24, 2011

### I like Serena

All right.
Here's the code for the system:
dY/dx = y
dy/dx = fdy(x, y) / fdx(x, y, t)
dt/dx = 1 / fdx(x, y, t)​

Code (Text):

private const double H = 0.01; // This will now be the x step instead of the t step

public double t = 0;
public double x = 0.81;
public double y = 0.055;
public double Y = 0;
// ;
public ArrayList tList;
public ArrayList xList;
public ArrayList yList;
public ArrayList YList;

public void Start()
{
tList = new ArrayList();
xList = new ArrayList();
yList = new ArrayList();
YList = new ArrayList();

for (double krok = 0; krok <= 500; krok++)
{
rungeStep(H, x, y, Y, t);
}
}

public void rungeStep(double H, double x, double y, double Y, double t)
{
double K1Y = H * y;
double K1y = H * dy(x, y) / fdx(x, y, t)
double K1t = H * 1 / dx(x, y, t);

double K2Y = H * (y + 0.5 * K1y);
double K2y = H * dy(x + 0.5 * H, y + 0.5 * K1y) / dx(x + 0.5 * H, y + 0.5 * K1y, t + 0.5 * K1t);
double K2t = H * 1 / dx(x + 0.5 * H, y + 0.5 * K1y, t + 0.5 * K1t);

double K3Y = H * (y + 0.5 * K2y);
double K3y = H * dy(x + 0.5 * H, y + 0.5 * K2y) / dx(x + 0.5 * H, y + 0.5 * K2y, t + 0.5 * K2t);
double K3t = H * 1 / dx(x + 0.5 * H, y + 0.5 * K2y, t + 0.5 * K2t);

double K4Y = H * (y + K3y);
double K4y = H * dy(x + H, y + K3y) / dx(x + H, y + K3y, t + K3t);
double K4t = H * 1 / dx(x + H, y + K3y, t + K3t);

double stepY = (K1Y + 2 * K2Y + 2 * K3Y + K4Y) / 6;
double stepy = (K1y + 2 * K2y + 2 * K3y + K4y) / 6;
double stept = (K1t + 2 * K2t + 2 * K3t + K4t) / 6;

this.x = x + H;
this.Y = Y + stepY;
this.y = y + stepy;
this.t = t + stept;
}

Last edited: Jul 24, 2011
20. Jul 25, 2011

### I like Serena

I just realized there's another differential equation (DE) that you can use.
The advantage is that you get exactly what you got before, but you also get the area under curve as a function of time (instead of x).

The DE is:
dx/dt = fdx(x, y, t)
dy/dt = fdy(x, y)
dY/dt = y fdx(x, y, t)​

That is, just one additional evaluation to keep track of the area under curve.
It will work without any extra conditions.