'********************************************************** '* PROGRAM CIRCUIT * '* Use Runge-Kutta method to solve a LRC circuit or equi- * '* valent damped mass-spring problem. * '* ------------------------------------------------------ * '* SAMPLE RUN: * '* * '* (Solve mass-spring system with: * '* F=75 N, M=50 kg, K=100 N (C=0.01), D=0.05) * '* * '* Specify [o]scillating or [n]o oscillating term... * '* o * '* Give V or F, L or M, R or D, C or 1/K (4L/C-R^2 > 0): * '* 75 50 0.05 0.01 * '* one period at t= 4.442879 s * '* Give t* 10 10 * '* * '* time q or d i or speed analytic q or d * '* 1.000000 0.7504 1.4993 0.6331 * '* 2.000000 2.2485 -0.0017 1.4628 * '* 3.000000 0.7483 -2.9955 1.0890 * '* 4.000000 -2.2440 0.0055 0.1436 * '* 5.000000 0.7555 5.9850 0.2223 * '* 6.000000 6.7320 -0.0150 1.1913 * '* 7.000000 0.7351 -11.9580 1.4147 * '* 8.000000 -11.2021 0.0378 0.5163 * '* 9.000000 0.7878 23.8922 0.0131 * '* 10.000000 24.6302 -0.0915 0.7537 * '* * '* (For a better accuracy, increase number of points). * '* * '* Basic version by J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '********************************************************** 'Program Circuit; DEFDBL A-H, O-Z DEFINT I-N 'main program '********************************************************* '* Use Runge-Kutta method to solve LRC circuit problems. * '* Variable equivalences with mass-and-spring problem: * '* V => force F * '* q => displacement x * '* i => velocity v * '* L => mass m * '* R => damping constant D * '* 1/C => spring constant K * '********************************************************} ' xi,q,XL,C,R,v,t,dt,t,tfinal: REAL ' xk1,xk2,xk3,xk4: REAL ' j,n: INTEGER; ' choice\$: CHAR ' Choose circuit type... CLS pi = 3.1415926535# PRINT INPUT " Specify [o]scillating or [n]o oscillating term... ", choice\$ IF choice\$ = "o" OR choice\$ = "O" THEN ' Ld^2q/dt^2+Rdq/dt+q/C=V 10 PRINT " Give V or F, L or M, R or 1/K, C or D (4L/C-R^2 > 0):" INPUT " ", vv, XL, R, C F = vv: XM = XL: XK = 1# / C: D = R PRINT " one period at t="; 4# * pi * XL / SQR(4# * XL / C - R * R); " s" PRINT " Give t_final and number of points:" INPUT " ", tfinal, n IF 4# * XL / C - R * R < 0# THEN GOTO 10 x = 0#: v = C * vv * R / 2# / XL: t = 0# 'Initial values dt = tfinal / n PRINT PRINT " time q or d i or speed analytic q or d" F1\$ = " ## ####.###### ####.###### ####.######" FOR j = 1 TO n GOSUB 1000 'MassAndSpring(x,v,R,1./C,L,vv,dt) t = t + dt PRINT USING F1\$; t; x; v; C * vv * (1# - EXP(-R * t / 2# / XL) * COS(SQR(4# * XL / C - R * R) * t / 2# / XL)) NEXT j ELSEIF choice\$ = "n" OR choice\$ = "N" THEN ' Ldi/dt+Ri=V (no oscillating term)... PRINT " Give V or F, L or M, R or 1/K:" INPUT " ", v, XL, R PRINT PRINT " time constant at t="; XL / R; " s" PRINT " Give t_final and number of points:" INPUT " ", tfinal, n xi = 0#: t = 0# 'Initial values dt = tfinal / n PRINT PRINT " time i or speed analytic i" F\$ = " ## ####.###### ####.######" FOR j = 1 TO n ' Runge-Kutta coefficients... xk1 = (v - R * xi) / XL xk2 = (v - R * (xi + xk1 * dt / 2#)) / XL xk3 = (v - R * (xi + xk2 * dt / 2#)) / XL xk4 = (v - R * (xi + xk3 * dt)) / XL ' Propagate solution... xi = xi + (xk1 + 2# * xk2 + 2# * xk3 + xk4) * dt / 6# t = t + dt PRINT USING F\$; t; xi; v / R * (1# - EXP(-R * t / XL)) NEXT j ELSE PRINT " No such choice. Try again..." END IF END '------------------------------------------------------- 500 'FUNCTION AofT(x,v,D,K,M,F:REAL):REAL; ' Calculate acceleration for "mass and spring" problem. AofT = -(-F + D * v1 + XK * x1) / XM RETURN '------------------------------------------------------------------------- 1000 'PROCEDURE MassAndSpring(VAR x,v:REAL; D,K,M,F,dt:REAL); ' Calculate motion for mass and spring problem with constant force term. ' xk1x,xk2x,xk3x,xk4x,xk1v,xk2v,xk3v,xk4v: REAL ' Runge-Kutta coefficients... xk1x = v x1 = x: v1 = v: GOSUB 500 xk2x = v + AofT * dt / 2# xk3x = v + AofT * dt / 2# xk4x = v + AofT * dt xk1v = AofT x1 = x + v * dt / 2#: v1 = v + xk1v * dt / 2# GOSUB 500: xk2v = AofT v1 = v + xk2v * dt / 2# GOSUB 500: xk3v = AofT x1 = x + v * dt: v1 = v + xk3v * dt GOSUB 500: xk4v = AofT ' Propagate solution... x = x + (xk1x + 2# * xk2x + 2# * xk3x + xk4x) * dt / 6# v = v + (xk1v + 2# * xk2v + 2# * xk3v + xk4v) * dt / 6# RETURN 'end of file circuit.bas