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


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

F = x

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

G = a*x + b*y

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
! 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"'
res = SYSTEMQQ("start " // PlotScriptFile)

end program RungeKutta

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.

  • $\begingroup$ besides, why did you add a cutoff on $V$? $\endgroup$ – Benoit Gaudeul Jul 10 '18 at 15:48

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.

  • $\begingroup$ Thanks, the problem was time steps. I should select very small time steps, at least $10^{-3}$. $\endgroup$ – Masoud Jul 11 '18 at 7:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.