Edwards-Penney Chapter 4, sections 4.1, 4.2, 4.3, 4.4 Topics, Definitions, Theorems, Examples 4.1 First order systems and applications ==== EXAMPLE. Forced system of undamped oscillators m1 x'' = - k1 x + k2 (y-x), m2 y'' = -k2(y-x) + f(t) INSTANCE m1=2, m2=1, k1=4, k2=2, f(t)=40sin(3t) 2 x'' = - 6x + 2y, y'' = 2x - 2y + 40 sin(3t) EXAMPLE. Recirculating brine tanks 20 x' = -6x + y, 20 y' = 6x - 3y x(t)=pounds of salt in tank 1 (100 gal) y(t)=pounds of salt in tank 2 (200 gal) x(0), y(0) = initial salt amounts in each tank t=minutes 20=inflow rate=outflow rate 0=inflow salt concentration FIRST ORDER SYSTEMS EXAMPLE. Convert x''' + 3x'' + 2x' - 5x = sin(2t) into a first order system using the position-velocity substitution x1=x, x2=x', x3=x''. ANSWER: x1' = x2, x2' = x3, x3' = 5x1 - 2x2 - 3x3 + sin(2t) EXAMPLE. Transform to a first order system x'' = x^3 + (x')^3 ANSWER: x1=x, x2=x'', x1'=x2, x2'=(x1)^3+(x2)^3 EXAMPLE. Transform to a first order system 2x'' = -6x + 2y, y'' = 2x - 2y + 40 sin(3t) ANSWER: x1=x,x2=x',x3=y,x4=y' ==> x1' = x2, x2' = -3x1 + x3, [a division by 2 needed] x3' = x4, x4' = 2x1 - 2x3 + 40 sin(3t) DYNAMICAL SYSTEMS The typical second order equation from circuit theory and mechanical oscillations is x'' + px' + qx = 0. It is realized as a 2-dimensional system x' = y, y' = -qx - py EXAMPLE. Solve x'=-2y, y'=x/2. ANSWER: x(t)=A cos(t) + B sin(t), y(t) = (-A/2) cos(t) + (B/2) sin(t). METHOD: Differentiate the first equation, substitute the second to get x''+x=0. Then x(t) is a linear combination of cos(t), sin(t). From the first equation, y=(-1/2)x', which implies the equation for y(t) by differentiation. This method is called the CAYLEY-HAMILTON SHORTCUT. THEOREM. [Cayley-Hamilton-Ziebur] A linear dynamical system u'=Au with 2x2 constant matrix A can be solved by finding the roots of the characteristic equation |A-rI|=0. Then x(t) and y(t) are linear combinations of the solution atoms from Euler's Theorem corresponding to the roots. EXAMPLE. Solve x'=-2y, y'=x/2, x(0)=2, y(0)=0. ANSWER: x(t)=2 cos(t), y(t)=sin(t) METHOD: Find A,B in the general solution from the two initial conditions at t=0. EXAMPLE. Solve x'=y, y'=2x+y ANSWER: x(t)=A exp(-t) + B exp(2t), y(t)=-a exp(-t)+2B exp(2t) METHOD: Cayley-Hamilton shortcut; A:=Matrix([[0,1],[2,1]]). Char equation r^2-r-2=0, roots r=-1,2, solution atoms are exp(-t), exp(2t). Then x(t)=linear combination. To find y, use the first equation, y=x'. EXAMPLE. Solve x'=-y, y'=1.01 x - 0.2 y, x(0)=0, y(0)=1. ANSWER: x(t)=exp(-t/10) sin(t), y(t)=(1/10)exp(-t/10)(sin(t)+10 cos(t)) METHOD: Cayley-Hamilton-Ziebur shortcut; A:=Matrix([[0,-1],[1.01,-0.2]]) with char equation r^2+0.2r+1.01=0 having complex roots 0.1+i, 0.1-i and corresponding Euler solution atoms exp(-t/10)cos(t), exp(-t/10)sin(t). Then x(t) is a linear combination of the solution atoms and y=-x' from the first equation. To find A,B in the general solution, use the two initial conditions at t=0. THEOREM 1. [Picard's Existence-Uniqueness for linear systems] Consider the system of differential equations u'=P(t)u+F(t) where u=, P(t) is nxn with continuous coefficients, F(t)= has continuous entries. Let b= be a constant vector and t=a a point in the interval J of continuity. Then the problem u' = P(t)u + F(t), u(a)=b, has a unique solution u(t) defined on J. 4.2 The Method of Elimination [Cayley-Hamilton-Ziebur Method] ==== EXAMPLE 1. Solve x' = 4x - 3y, y' = 6x - 7y, x(0)=2, y(0)=-1. SOLUTION: Elimination is equivalent to using Cayley-Hamilton-Ziebur. This example shows why. The characteristic equation is identical to the operational determinant, with roots r=2,-5 and Euler solution atoms e^(2t), e^(-5t)$. Then we have the conclusion of C-H_Z, that x(t) and y(t) are linear combinations of the Euler atoms: x(t) = a1 e^(2t) + a2 e^(-5t), y(t) = b1 e^(2t) + b2 e^(-5t). Exactly as in C-H-Z, these equations are substituted into the first differential equation, then coefficients of atoms are matched to get equations for a1,a2,b1,b2. ANSWER: x(t) = (3/2) b1 e^(2t) + (1/3) b2 e^(-5t), y(t) = b1 e^(2t) + b2 e^(-5t). EXAMPLE 2. Solve x'' + 3x - y = 0, -2x + y'' + 2y = 0. SOLUTION: The system can be represented as a 4x4 linear system of the form Au=0. The operational determinant is identical to |A - rI|, which is the characteristic equation of matrix A. According to Cayley-Hamilton-Ziebur, the components x,y are linear combinations of the Euler atoms found from the four roots i, -i, 2i, -2i of |A-rI|=0. The shortcut applies again, using the equations x(t) = a1 cos(t) + a2 sin(t) + a3 cos(2t) + a4 sin(2t), y(t) = b1 cos(t) + b2 sin(t) + b3 cos(2t) + b4 sin(2t). The plan is to substitute these equations into the first differential equation, then match coefficients of atoms to obtain formulas for b1 to b4 in terms of a1 to a4. ANSWER: x(t) = a1 cos(t) + a2 sin(t) + a3 cos(2t) + a4 sin(2t), y(t) = 2a1 cos(t) + 2a2 sin(t) - b1 cos(2t) - b2 sin(2t). DEFINITION: The NATURAL FREQUENCIES of oscillation are the trig frequencies 1, 2 in the example. See Example 1, section 4.1. 4.3 Numerical Methods for Systems ==== DE: du/dt = F(t,u), u=vector, t=scalar EULER METHOD. The vector formulas for iteration are t[n+1] = t[n] + h u[n+1] = u[n] + h F(t[n],u[n]) Given are u[0], t[0], h. The symbol F is the vector function defining the differential equation du/dt = F(t,u). IMPROVED EULER METHOD. HEUN METHOD. The vector formula for iteration is t[n+1] = t[n] + h Time variable y[n+1] = u[n] + h * F(t[n],u[n]) Predictor u[n+1] = u[n] + h/2 * F(t[n],u[n]) + h/2 * F(t[n+1],y[n+1]) Corrector Given are u[0], t[0], h. The symbol F is the vector function defining the differential equation du/dt = F(t,u). Intermediate symbol y[n] is used only to pass the predictor value to the corrector. RK4 METHOD. The vector formula for iteration has five lines: t[n+1] = t[n] + h k1 = F(t[n],u[n]) Predictor 1 k2 = F(t[n]+h/2,u[n] + k1*h/2) Predictor 2 k3 = F(t[n]+h/2,u[n] + k2*h/2) Predictor 3 k4 = F(t[n+1],u[n] + k3*h) Predictor 4 u[n+1] = u[n] + h(k1 + 2*k2 + 2*k3 + k4)/6 Corrector Given are u[0], t[0], h. The symbol F is the vector function defining the differential equation du/dt = F(t,u). Intermediate predictors k1, k2, k3, k4 are only used for the corrector, the fifth formula, their values forgotten after one use. EXAMPLE 1. Solve numerically the system du/dt = F(t,u), u[0] = (3,6) where u=(x(t),y(t)), F(t,u) = (3x-2y,5x-4y). Compare the calculations to the exact solution x(t)=2 e^(-2t) + e^t, y(t) = 5e^(-2t) + e^t using the Euler Method to find values for u(0.1), u(0.2) and u(0.3). ANSWERS: u[0]=(3,6), t[0]=0, h=0.1. A hand calculator is useful. Then u(0.1) approx u[1], u(0.2) approx u[2] and u(0.3) approx u[3]. u[1] = u[0] + h*F(0,u[0]) = (2.7,5.1) u[2] = u[1] + h*F(h,u[1]) = (2.4,4.2) u[3] = u[2] + h*F(2h,u[2]) = (2.58,4.62) EXAMPLE 2. Solve numerically the system du/dt = F(t,u), u[0] = (0,1) where u=(x(t),y(t)), F(t,u) = (y,-x). Compare the calculations to the exact solution x(t)=sin(t), y(t) = cos(t). Use the vector RK4 Method to find values for u(t) from t=0 to t=5 in steps h=0.5. ANSWER: See table 4.3.1 in the textbook. EXAMPLE 3. Solve numerically the Lunar Lander system du/dt = F(t,u), u[0] = (1781870,-450) where u=(x(t),y(t)), F(t,u) = (y,4 - 4.9044*10^12 /x^2). Use the RK4 Method with h=0.1. COMPUTER CODE: See the section on comets and spacecraft. Maple code for an answer check: des:=diff(x(t),t)=y(t),diff(y(t),t)=4 - 4.9044*10^12 /x(t)^2; ic:=x(0)=1781870,y(0)=-450; p:=dsolve({des,ic},{x(t),y(t)},numeric); p(100); p(200); EXAMPLE 4. Baseball problem. The equations are du/dt = F(t,u) u(0)=u[0]=(0,0,80*sqrt(3),80) h=0.1, 0.05 and 0.025 for multiple computations F(t,u) = (f1,f2,f3,f4) f1 = u3 f2 = u4 f3 = -c*u3*sqrt(u3^2+u4^2) f4 = -g - c*u4*sqrt(u3^2+u4^2) c = 0 No Air Resistance c = 0.0025 Nonlinear Newton drag g = 32 ANSWERS: For c=0 the exact solution is x(t)=80*sqrt(3)*t, y(t)=80t-16*t^2. Here's a quick answer check: des0:=diff(x(t),t,t)=0,diff(y(t),t,t)=-32; ic:=x(0)=0,y(0)=0,D(x)(0)=80*sqrt(3),D(y)(0)=80; dsolve({des0,ic},{x(t),y(t)}); For c=0.0025, no exact solution was attempted. Instead, we code for maple: c=0.0025; des1:=diff(x(t),t,t)=-c*sqrt((D(x)(t))^2+(D(y)(t))^2)*D(x)(t), diff(y(t),t,t)=-32-c*sqrt((D(x)(t))^2+(D(y)(t))^2)*D(y)(t); ic:=x(0)=0,y(0)=0,D(x)(0)=80*sqrt(3),D(y)(0)=80; p:=dsolve({des1,ic},{x(t),y(t)},numeric); An answer check for Table 4.3.5 columns 1,2 can be made by printing p(t), for t=0 to t=4 step 0.5, e.g., p(4.0) is [t = 4., x(t) = 339.671376642473, diff(x(t), t) = 54.7872713137340, y(t) = .908591100591405, diff(y(t), t) = -56.4425175342350] The velocity is given by the formula V:=t->sqrt(rhs(p(4)[3])^2+rhs(p(4)[5])^2); with V(4) equal to 78.66004629802305, agreeing with Table 4.3.5. ===== end