fourth order Runge–Kutta in C#  two differential equationsby Itosu Tags: differential, equations, order, runge–kutta 

#1
Jul2011, 06:07 PM

P: 19

Hello!
I am trying to solve two differential, nonlinear equations in C# using fourthorder RungeKutha method. I implemented them in C#:
Do you have any ideas? Thank you in advance! 



#2
Jul2011, 11:08 PM

HW Helper
P: 6,189

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. Anyway, your code




#3
Jul2111, 06:51 AM

P: 19

Newer version below:




#4
Jul2111, 07:46 AM

P: 19

fourth order Runge–Kutta in C#  two differential equations
The newest version is like this:
any comments? what about "t"? 



#5
Jul2111, 09:24 AM

HW Helper
P: 6,189

It should be like this:




#6
Jul2111, 09:35 AM

P: 19

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
Jul2111, 10:28 AM

HW Helper
P: 6,189

It should be:
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. 



#8
Jul2111, 12:21 PM

P: 19

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



#9
Jul2111, 01:00 PM

HW Helper
P: 6,189

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(h^{5}). 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
Jul2211, 08:34 AM

P: 19

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? and.. thank you in advance! 



#11
Jul2311, 09:55 AM

HW Helper
P: 6,189

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%...Cotes_formulas 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 f_{0}, f_{1}, f_{2}, ..., f_{n} are your function values with a step size H between them. You have for instance the options: Trapezoid rule (error is O(H^{3}) which is less accurate than RungeKutta). n * H * (f_{0}/2 + f_{1} + f_{2} + f_{3} + ... + f_{(n1)} + f_{n}/2) Simpson's rule with n even (error is O(H^{5}) which is more accurate than RungeKutta): (n * H / 6) * (f_{0} + 4 f_{1} + 2 f_{2} + 4 f_{3} + 2 f_{4} + ... + 2 f_{(n1)} + f_{n}) 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. 



#12
Jul2311, 11:34 AM

Engineering
Sci Advisor
HW Helper
Thanks
P: 6,340

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
Jul2411, 05:48 AM

P: 19

What do you mean? Can you explain?




#14
Jul2411, 09:44 AM

HW Helper
P: 6,189

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 RungeKutta 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
Jul2411, 10:21 AM

P: 19

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
Jul2411, 11:58 AM

HW Helper
P: 6,189

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)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 (x_{0}, y_{0}, 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=x_{0} 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 = yDifferentiation 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 = yThere 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
Jul2411, 12:13 PM

P: 19

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


Register to reply 
Related Discussions  
Solving Second Order Differential Equations using Runge Kutta  Calculus & Beyond Homework  7  
2nd Order RungeKutta: 2nd Order Coupled Differential Equations  Calculus & Beyond Homework  3  
Fourthorder Runge–Kutta method? =(  Calculus & Beyond Homework  2  
Fourth Order RungeKutta Method  Programming & Computer Science  10 