rk4: Classical Runge-Kutta

Description Usage Arguments Details Value Note References See Also Examples

Description

Classical Runge-Kutta of order 4.

Usage

1
2
3
rk4(f, a, b, y0, n)

rk4sys(f, a, b, y0, n)

Arguments

f

function in the differential equation y' = f(x, y);
defined as a function R \times R^m \rightarrow R^m, where m is the number of equations.

a, b

endpoints of the interval.

y0

starting values; for m equations y0 needs to be a vector of length m.

n

the number of steps from a to b.

Details

Classical Runge-Kutta of order 4 for (systems of) ordinary differential equations with fixed step size.

Value

List with components x for grid points between a and b and y an n-by-m matrix with solutions for variables in columns, i.e. each row contains one time stamp.

Note

This function serves demonstration purposes only.

References

Fausett, L. V. (2007). Applied Numerical Analysis Using Matlab. Second edition, Prentice Hall.

See Also

ode23, deval

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
##  Example1: ODE
# y' = y*(-2*x + 1/x) for x != 0, 1 if x = 0
# solution is x*exp(-x^2)
f <- function(x, y) {
	if (x != 0) dy <- y * (- 2*x + 1/x)
	else        dy <- rep(1, length(y))
	return(dy)
}
sol <- rk4(f, 0, 2, 0, 50)
## Not run: 
x <- seq(0, 2, length.out = 51)
plot(x, x*exp(-x^2), type = "l", col = "red")
points(sol$x, sol$y, pch = "*")
grid()
## End(Not run)

##  Example2: Chemical process
  f <- function(t, u) {
    u1 <- u[3] - 0.1 * (t+1) * u[1]
    u2 <- 0.1 * (t+1) * u[1] - 2 * u[2]
    u3 <- 2 * u[2] - u[3]
    return(c(u1, u2, u3))
  }
u0 <- c(0.8696, 0.0435, 0.0870)
a <- 0; b <- 40
n <- 40
sol <- rk4sys(f, a, b, u0, n)
## Not run: 
matplot(sol$x, sol$y, type = "l", lty = 1, lwd = c(2, 1, 1),
        col = c("darkred", "darkblue", "darkgreen"),
        xlab = "Time [min]", ylab= "Concentration [Prozent]",
        main = "Chemical composition")
grid()
## End(Not run)


Search within the pracma package
Search all R packages, documentation and source code

Questions? Problems? Suggestions? or email at ian@mutexlabs.com.

Please suggest features or report bugs with the GitHub issue tracker.

All documentation is copyright its authors; we didn't write any of that.