Learn Runge-Kutta for Solving Differential Equations

[saltydog]

I'd like to learn how to use Runge-Kutta to solve IDEs. My understanding is that this technique can be used for the following:

\frac{dy}{dt}=F(t,y)+\int_0^t G[t,s,y(s)]ds

So I start with the following simple version:

\frac{dy}{dt}=ay-by^2+K\int_0^t y(s)ds

This is a good start because this equation can be differentiated to yield:

y(y^{''})=(y^{'})^{2}-by^{'}y^{2}+Ky^3

with:

y^{'}(0)=ay(0)-b[y(0)]^2

And so I can numerically solve the resulting ODE to get a precise answer that I can check against the Runge-Kutta version. Well, my results aren't "close" enough. I've included two plots: The first plot is the Runge-Kutta (dark line) against the ODE solution. The second plot gives the absolute difference between the two. The third plot is the back-substitution of the numerical results of the ODE back into the IDE. As you can see, the values are no more than 10^{-5} which is close.

However, the difference between the ODE results and Runge-Kutta of the IDE is the same type even with smaller increments of t and different values of a and b and K.

Perhaps I'm implementing Runge-Kutta incorrectly. Can anyone give me suggestions?

Thanks,
Salty

 

[saltydog]

Ok, I'm making progress: Instead of using Runge-Kutta, I've dropped back to just trying an Euler approximation first (still wish a Runge-Kutta).

To the value of the slope calculated by the Euler method, I just add the running value of the developing integral (calculated at each increment of t). This gives me a better approximation which improves with smaller increment in t. Using the function calculated by the Euler method, I determine the LHS and RHS of the IDE and take the absolute value of the difference. The attached plot shows one such run for delta t set at 0.001. As you can see, the difference is not greater than 0.001.

I'm not doing something right with Runge-Kutta else I'd expect a better match to the actual result.

 

[saltydog]

I made a change to the Euler method I'm using to evaluate the IDE. I've attached a comparison plot of the solution calculated by the ODE and a solution calculated directly from the IDE using a modified Euler method. As you can see, the plots overlap nicely but the Euler method is still not precise enough as seen by the back substitution into the IDE. The second plot shows the difference: LHS-RHS of the IDE. Some values are as large as 10^{-3} with delta t set at 10^{-3}. A slight improvement is obtained with smaller delta t.

Really I'd like to perfect a very accurate numerical approach as I'd like to study some coupled equations over long time periods. The on-line resource is lacking. If no one has an idea, I'll just jump in my truck next week and try the old-fashioned approach: go to the library.

 

[saltydog]

I think I may have a working Runge-Kutta method to solve this IDE. Using the technique I designed for the Euler method, I carried over to the Runge-Kutta method. Really, the Runge-Kutta is a version of Euler in which a weighted average of 4 slopes are calculated as opposed to one. Using the same equation with step size of 0.001 in the interval of (0,20), I was able to lower the difference of the LHS-RHS of the IDE to 10^{-7}. The attached plot shows this difference.

Really, I need to experment with different IDEs to see how the error varies and also still need to compare my method with whatever published methods are available. It's a nice ride to the campus.

 

[saltydog]

I've been working with:

y^{'}=ay-by^{2}-y(1+t)\int_0^ty(s)ds

Interesting that this one can also be differentiated yielding the very impressive looking:

y^2y^{'''}-y^{''}(3yy^{'}-by^{3})+2(y^{'})^{3}=-(1+t)y^{3} y^{'}-2y^{4}

Mathematica solves this easily and again I can compare this with my numerical results for the IDE.

I've only looked at it using an Euler method with ok results (within 10^{-3} for 0.01 step size.

I'll convert it to Runge-Kutta and also begin looking at IDEs which can't be so nicely differentiated such as those with kernels that can't be eliminated by differentiation. That (and systems) is really why I want to have a good numerical method.