## Re: [Maxima-discuss] How to expand physical constants?

 Re: [Maxima-discuss] How to expand physical constants? From: Robert Dodier - 2014-06-04 22:05:14 ```On 2014-06-01, Gunter Königsmann wrote: > One example would be a damped oscillator: That's an interesting problem. Although it's already been solved, I feel obliged to throw my hat in the ring, as they say. I tried a couple of different approaches, and as it turns out, it worked best to solve the equation and substitute dimensional values into the result. So that's what I'll show here. (%i1) display2d : false \$ (%i2) load (ezunits) \$ (%i3) Values : [C = 1e-9 ` farad, L = 1e-6 ` henry, R = 100.0 ` ohm] \$ (%i4) g1: U_R(t)=R*I(t) \$ (%i5) g2: U_L(t)=L*'diff(I(t),t) \$ (%i6) g3: I(t)=C*'diff(U_C(t),t) \$ (%i7) g4: U_C(t)+U_R(t)+U_L(t)=0 \$ (%i8) dgl:ev(g4,[g1,g2,g3]) \$ I wonder if U_R, U_L, and U_C, as defined, all have the same dimensions. (%i9) [dimensions (ohm * A), dimensions (henry * A / s), dimensions (s * A / farad)]; (%o9) [length^2*mass/(current*time^3),length^2*mass/(current*time^3), length^2*mass/(current*time^3)] OK, that's reassuring. Now considering the asksign that we're going to get from desolve, (%i10) C*(C*R^2 - 4*L), Values; (%o10) (1.0E-9 ` farad)*(1.0E-5 ` farad*ohm^2 +(-3.9999999999999997E-6) ` henry) (%i11) foo : expand (%); (%o11) 1.0E-14 ` farad^2*ohm^2+(-4.0E-15) ` farad*henry Are those the same dimensions? I hope so. (%i12) dimensions (%); (%o12) 2*time^2 OK, reassuring again. Convert to units for time^2. (%i13) foo `` s^2; ezunits: computing conversions to base units; may take a moment. (%o13) 6.000000000000001E-15 ` s^2 So that's positive. Now solve the differential equation. (%i14) foo : desolve (ev (dgl, nouns), U_C(t)); Is C*(C*R^2-4*L) positive, negative or zero? p; (%o14) U_C(t) = %e^-(t*R/(2*L))*((2*C*L *(U_C(0)*C*R +('at('diff(U_C(t),t,1),t = 0))*C*L) -U_C(0)*C^2*L*R) *sinh(t*sqrt(C*(C*R^2-4*L))/(2*C*L)) /sqrt(C*(C*R^2-4*L)) +U_C(0)*C*L *cosh(t*sqrt(C*(C*R^2-4*L))/(2*C*L))) /(C*L) Subsitute in initial values. Specify U_C(0) in appropriate units. (%i15) atvalue ('diff(U_C(t), t, 1), t = 0, 0) \$ (%i16) atvalue (U_C (t), t = 0, 10 ` ohm * A) \$ (%i17) rhs (foo), at; (%o17) %e^-(t*R/(2*L))*(10*C^2*L*R*sinh(t*sqrt(C*(C*R^2-4*L))/(2*C*L)) /sqrt(C*(C*R^2-4*L)) +10*C*L*cosh(t*sqrt(C*(C*R^2-4*L))/(2*C*L))) /(C*L) ` ohm*A Substitute parameters C, L, and R. (%i18) bar : subst (Values, %); (%o18) ((9.999999999999999E+14 ` 1/henry)*%e^((-5.0E+7*t) ` ohm/henry) ` 1 /farad) *((1.0E-21*sinh(sqrt((1.0E-9 ` farad)*(1.0E-5 ` farad*ohm^2 +(-3.9999999999999997E-6) ` henry)) *(4.999999999999999E+14*t ` 1/(farad*henry))) ` farad^2*henry*ohm) /sqrt((1.0E-9 ` farad)*(1.0E-5 ` farad*ohm^2 +(-3.9999999999999997E-6) ` henry)) +1.0E-14*cosh(sqrt((1.0E-9 ` farad)*(1.0E-5 ` farad*ohm^2 +(-3.9999999999999997E-6) ` henry)) *(4.999999999999999E+14*t ` 1/(farad*henry))) ` farad*henry) ` ohm*A Hmm ... there are lots of combinations of units, which are (I hope) all equivalent. Let's change those to the equivalent basic units. (%i19) matchdeclare ([a%, b%], all) \$ (%i20) defrule (r1, a% ` b%, a% ` fundamental_units (b%)) \$ (erm ... certain amount of flailing around here ... OK, back on track.) (%i31) applyb1 (bar, r1); (%o31) (9.999999999999999E+14 ` s^2*A^2/(kg*m^2)) *%e^((-5.0E+7*t) ` 1/s) *(1.290994448735805E-14*sinh(3.872983346207417E+7*t ` 1/s) +1.0E-14*cosh(3.872983346207417E+7*t ` 1/s)) ` kg^2*m^4/(s^5*A^3) Time is measured in seconds. Let's represent time as a dimensionless variable t1. (%i32) %, t = t1 ` s; (%o32) %e^((-5.0E+7*t1) ` 1)*(9.999999999999999E+14 *(1.290994448735805E-14 *sinh(3.872983346207417E+7*t1 ` 1) +1.0E-14*cosh(3.872983346207417E+7*t1 ` 1)) ` s^2*A^2/(kg*m^2)) ` kg^2*m^4/(s^5*A^3) Encourage simplification of (stuff) ` 1 to (stuff). (%i33) ''%; (%o33) 9.999999999999999E+14*%e^-(5.0E+7*t1) *(1.290994448735805E-14 *sinh(3.872983346207417E+7*t1) +1.0E-14*cosh(3.872983346207417E+7*t1)) ` kg*m^2/(s^3*A) Make a nice plot. (%i34) plot2d (qty (%), [t1, 1e-10, 1e-6]) \$ OK, well, that wasn't so painful. Kind of fun, actually. best, Robert Dodier ```