/********************************************************* * 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_final and number of points: * * 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 * * * * Press any key to continue * * (For a better accuracy, increase number of points). * * * * ------------------------------------------------------ * * Reference: "Problem Solving with Fortran 90 By David * * R. Brooks, Springer-Verlag New York, 1997" * * * * C++ Version By J-P Moreau, Paris. * * (www.jpmoreau.fr) * *********************************************************/ #include #include //------------------------------------ double AofT(double x,double v,double D,double K, double M, double F) { // Calculate acceleration for "mass and spring" problem. return -(-F+D*v+K*x)/M; } //--------------------------------------------- void MassAndSpring(double *x,double *v, double D,double K, double M,double F,double dt) { // Calculate motion for mass and spring problem with constant force term. double k1_x,k2_x,k3_x,k4_x,k1_v,k2_v,k3_v,k4_v; // Runge-Kutta coefficients... k1_x=*v; k2_x=*v+AofT(*x,*v,D,K,M,F)*dt/2.0; k3_x=*v+AofT(*x,*v,D,K,M,F)*dt/2.0; k4_x=*v+AofT(*x,*v,D,K,M,F)*dt; k1_v=AofT(*x,*v,D,K,M,F); k2_v=AofT(*x+(*v*dt/2.0),*v+k1_v*dt/2.0,D,K,M,F); k3_v=AofT(*x+(*v*dt/2.0),*v+k2_v*dt/2.0,D,K,M,F); k4_v=AofT(*x+(*v*dt),*v+k3_v*dt,D,K,M,F); // Propagate solution... *x=*x+(k1_x+2.0*k2_x+2.0*k3_x+k4_x)*dt/6.0; *v=*v+(k1_v+2.0*k2_v+2.0*k3_v+k4_v)*dt/6.0; } // MassAndSpring() void main() { /******************************************************** * 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 * ********************************************************/ double i,q,L,C,R,V,t,dt,t_final; double k1_i,k2_i,k3_i,k4_i,pi=3.1415926535; int j,n; char choice; // Choose circuit type... printf("\n Specify [o]scillating or [n]o oscillating term... "); scanf("%c",&choice); if (choice=='o' || choice=='O') { // Ld^2q/dt^2+Rdq/dt+q/C=V do { printf(" Give V or F, L or M, R or 1/K, C or D (4L/C-R^2 > 0):\n"); printf(" "); scanf("%lf %lf %lf %lf",&V,&L,&R,&C); printf(" one period at t=%f s\n", 4*pi*L/sqrt(4.0*L/C-R*R)); printf(" Give t_final and number of points:\n"); printf(" "); scanf("%lf %d", &t_final,&n); } while (4.0*L/C-R*R <= 0.0); q=0.0; i=C*V*R/2.0/L; t=0.0; //Initial values dt=t_final/n; printf("\n time q or d i or speed analytic q or d\n"); for (j=1; j<=n; j++) { MassAndSpring(&q,&i,R,1.0/C,L,V,dt); t=t+dt; printf(" %12.6f%12.4f%12.4f%12.4f\n", t, q, i, C*V*(1.0-exp(-R*t/2.0/L) *cos(sqrt(4.0*L/C-R*R)*t/2.0/L))); } } else if (choice=='n' || choice=='N') { // Ldi/dt+Ri=V (no oscillating term)... printf(" Give V or F, L or M, R or 1/K:\n"); printf(" "); scanf("%lf %lf %lf", &V,&L,&R); printf(" time constant at t=%f s\n", L/R); printf(" Give t_final and number of points:\n"); printf(" "); scanf("%lf %d", &t_final, &n); i=0.0; t=0.0; //Initial values dt=t_final/n; printf("\n time i or speed analytic i\n"); for (j=1; j<=n; j++) { // Runge-Kutta coefficients... k1_i=(V-R*i )/L; k2_i=(V-R*(i+k1_i*dt/2.0))/L; k3_i=(V-R*(i+k2_i*dt/2.0))/L; k4_i=(V-R*(i+k3_i*dt) )/L; // Propagate solution... i=i+(k1_i+2.0*k2_i+2.0*k3_i+k4_i)*dt/6.0; t=t+dt; printf(" %12.6f,%12.4f%12.4f\n", t, i, V/R*(1.0-exp(-R*t/L))); } } else printf(" No such choice. Try again...\n"); printf("\n"); } // end of file circuit.cpp