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

## Main Question or Discussion Point

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:
 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?

#### Attachments

• 3.8 KB Views: 808

Related Differential Equations News on Phys.org
I like Serena
Homework Helper
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:
        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:
        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);

Last edited:
The newest version is like this:
Code:
        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:
I like Serena
Homework Helper
The newest version is like this:

and looks like the one in Simulink :)
Yes, I have some comments. :)
It should be like this:

Code:
        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!

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 :)

I like Serena
Homework Helper
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 :)
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:
        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:
You are totally correct. Thank you!
Do you know of any possibilities to calculate an error of such solution?

I like Serena
Homework Helper
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.

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?

I like Serena
Homework Helper
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:
AlephZero
Homework Helper
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?
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?

What do you mean? Can you explain?

I like Serena
Homework Helper
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.

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!

I like Serena
Homework Helper
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?

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?
Yes.

Can I assume that x strictly increases with time?
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 :)

#### Attachments

• 5.8 KB Views: 686
I like Serena
Homework Helper
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?

I like Serena
Homework Helper
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:
        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:
I like Serena
Homework Helper
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.

Could you explain in more details where dY/dt = y fdx(x, y, t) came from?

And to calculate error do I understand it correctly that I have to run the same programme with H=0.01/2 and compare results?

Thank you one more time!

I like Serena
Homework Helper

Could you explain in more details where dY/dt = y fdx(x, y, t) came from?

And to calculate error do I understand it correctly that I have to run the same programme with H=0.01/2 and compare results?

Thank you one more time!
The chain rule says that:
d/dt Y(x(t)) = dY/dx dx/dt​

And since Y is the integral of y with respect to x, y is the derivative of Y with respect to x, so:
dY/dx = y​

Combining the two gives:
dY/dt = y dx/dt​

Perhaps more intuitive is thinking of rectangles with an infinitesimal width dx.
When x increases to x+dx, the area under curve increases with dY = y dx, which is the surface of an infinitesimal rectangle.
This happens between time t and t+dt.
So dY/dt = y dx/dt.

And yes, running the program again with half of the step size will give you an estimation of the error.
Of course, when you have done that, this result will be the more accurate one. :)

Last edited:
Thank you. I understand everything (are you a teacher? :) but one thing:
And since Y is the integral of y with respect to x,
How do we know that? Why with respect to x, not to t?

And as far as error is concerned - I start with half smaller step size, and then substract results?

Thank you again!

I like Serena
Homework Helper
Thank you. I understand everything (are you a teacher? :) but one thing:

How do we know that? Why with respect to x, not to t?

And as far as error is concerned - I start with half smaller step size, and then substract results?

Thank you again!
Professionally I'm not a teacher but an engineer, but I do give private lessons. :)
And as you can see, I've become a homework helper on PF.

You asked for the area under curve of y with respect to x.
In a previous post I have defined Y to be this area under curve (with respect to x).
Isn't that what you wanted?

As for the error.
When you calculate the result twice with the 2 different step sizes, you'll find in each case a final (t, x, y, Y).
You can interpret your error in different ways.
Since you mentioned that the area under curve is carried forward, I would interpret the error as the different areas under curve at the final time t.

Btw, I'm still interested in what you are measuring.
What do your symbols x, y and Y represent?

I am full of admiration that you find a time to support us on the forum :)

x - glucose level
y - insuline level
Y - AUC of the glucose level function
I need AUC to calculate http://en.wikipedia.org/wiki/Glycemic_index" [Broken]

So, as I need AUC under x (as in previous image), on vertical axis there is x, on horizontal time, shouldn't we estimate the integral of x?

Last edited by a moderator: