solving RLC circuit equation using RK4

Question

I'm trying to solve RLC circuit differential equation numerically using 4th order Runge-Kutta method. $$LC\frac{\mathrm{d}^2V}{\mathrm{d}t^2}+RC\frac{\mathrm{d}V}{\mathrm{d}t}+V=0$$ First, I convert this equation to two first order equations with this change of variable: $$ U\equiv\frac{\mathrm{d}V}{\mathrm{d}t}\quad\Rightarrow\quad\frac{\mathrm{d}^2V}{\mathrm{d}t^2}=\frac{\mathrm{d}U}{\mathrm{d}t} $$ then $$ \begin{cases} \frac{\mathrm{d}V}{\mathrm{d}t}=U\equiv F(t,V,U) \\ \frac{\mathrm{d}U}{\mathrm{d}t}=-\dfrac{R}{L}U-\frac{1}{LC}V\equiv G(t,V,U) \end{cases} $$ I've written the following code in fortran90 to solve this problem. My problem is that, the code works correctly for capacities below about 0.00002, but when I input $C=0.000002$ the numbers grow hugely to the orders like $10^{28}$. I don't know why does this behaviour appear? I checked the code line by line many times, and I'm sure it is completely correct. The problem arises when I input very small $C$ or large $R$.

Could anyone please help me find out the problem?

This is my code in fortran90:

program RungeKutta
USE IFPORT
implicit none
! Constants
CHARACTER(50), PARAMETER :: outputFile = "output.dat", PlotConfigFile = "plotConfig.txt",  PlotScriptFile = "plotScript.bat"
real, parameter :: h = 1e-2, C = 2e-5, L = 0.5, R = 100.0, ti = 0.0, tf = 0.1
integer, parameter :: n = (tf - ti)/h
! Variables
real :: V(0:n), U0 = 0.0, U, a = -R/L, b = -1.0/(L*C), o6 = 1.0/6, Axis(0:n) = 0.0
real :: k0, k1, k2, k3, m0, m1, m2, m3
integer :: i

! Body of RungeKutta
V(0) = 1e2 ! initial condition

do i = 0, n-1
    k0 = h * F(U0)
    m0 = h * G(U0, V(i))

    k1 = h * F(U0 + m0/2)
    m1 = h * G(U0 + m0/2, V(i) + k0/2)

    k2 = h * F(U0 + m1/2)
    m2 = h * G(U0 + m1/2, V(i) + k1/2)

    k3 = h * F(U0 + m2)
    m3 = h * G(U0 + m2, V(i) + k2)

    V(i+1) = V(i) + o6*(k0 + 2*k1 + 2*k2 + k3)
    U = U0 + o6*(m0 + 2*m1 + 2*m2 + m3)

    if (abs(V(i+1)) < 1e-6) V(i+1) = 0.0

    Axis(i+1) = Axis(i) + h

    U0 = U
end do

call WriteOutputData()

CONTAINS

real FUNCTION F(x)
implicit none
real, intent(in) :: x

F = x
END FUNCTION

real FUNCTION G(x,y)
implicit none
real, intent(in) :: x, y

G = a*x + b*y
END FUNCTION

SUBROUTINE WriteOutputData()
implicit none
integer :: i
logical :: res
! Write voltages to output file
open(1, file = outputFile)
do i = 0, n
    write(1,'(F,1x,F)') Axis(i), V(i)
end do
close(1)
! Write plot config file
open(2, file = PlotConfigFile)
write(2,'(A)') 'set xlabel "Time (s)"'
write(2,'(A)') 'set ylabel "Voltage (V)"'
write(2,'(A)') 'plot "C:\\Users\\MASOUD\\Desktop\\numerical problems\\programs\\Runge-Kutta\\Runge-Kutta\\output.dat" with lines lt 3 lc 0 title "voltage"'
close(2)
res = SYSTEMQQ("start " // PlotScriptFile)
END SUBROUTINE

end program RungeKutta

Answer 1

For large $R/L$ or small $LC$ you system is stiff and you may quit the region of stability of RK4. In that case, you should turn to implicit methods or adapt your time step.

Answer 2

Your concrete equation with the data in the code is $$ V''+100V'+100000V=0 $$ The stability of RK4 was already mentioned in the other answer which would require step sizes smaller $2/\sqrt{LC}\sim$1e-3 for sensible results and about $10^{-3}/\sqrt{LC}\sim$1e-6 to 1e-5 for optimally correct results in 64bit double precision floating point.

One can also see this problem from the perspective of the sampling theorem. The first derivative term can be considered as a small perturbation, then the remaining equation is an oscillation with frequency $\sqrt{LC}=\sqrt{10^5}$. To get the wave form with any kind of reliability you will need (minimally 2 but better) 6 to 10 samples per period, which also gives the range of 2e-3 as an upper bound for the step size.