1Numerical Solutions of Systems of Numerical Solutions of Systems of Ordinary Differential EquationsOrdinary Differential EquationsLarry CarettoMechanical Engineering 501ABSeminar in Engineering AnalysisNovember 16, 20042Outline• Review last class– Stability of numerical solutions– Step size variation for error control– Multistep methods with constant and variable step size• Systems of equations• Definition of stiff systems• Algorithms for stiff systems• Finding if system is stiff3Systems of Equations• Any problem with one or more ODEs of any order can be reduced to a system for first order ODEs• In a previous lecture, we reduced a system of two second order equations to a system of four first order equations• We can reduce any system of equations into an equivalent system of first order equations• The sum of the orders is constant4Reduction of Order Example• Define variables y3= dy1/dt and y4= dy2/dt • Then dy3/dt = d2y1/dt2and dy4/dt = d2y2/dt2• Have four simultaneous first-order ODEs0022231222222121121212=++−=−++ ymkkymkdtydymkymkkdtyd),(),(),(),(422231224321211213242131yyyytfymkkymkdtdytfymkymkkdtdytfydtdytfydtdy=+−==++−=====5System of Equation Notation• Write system as vector equation),( yfyxdxd=• Individual components),,,,(),(21 NiiiyyyxfxfdxdyK== y• Each fimay depend on x (or t) and all yj),(),(),(),(422231224321211213242131yyyytfymkkymkdtdytfymkymkkdtdytfydtdytfydtdy=+−==++−=====6Review Example• Started with two second order equations (total order = 2 + 2 = 4)• Transformed into four first order equations (total order = 4 * 1 = 4)• Initial conditions on first derivatives of original variables in second order equation needed to get initial conditions on new variables like y3= dy1/dt• Look at numerical algorithms27Solving Simultaneous ODEs• Apply same algorithms used for single ODEs– Must apply each step and substep to all equations in system– Key is having consistent x and y values in determination of fi(x,y)– All yivalues in y must be available at the same x point when computing the fi– E.g., in Runge-Kutta we must evaluate k1for all equations before finding k28Runge-Kutta for ODE System– y(n)is vector of dependent variables at x = xn– k(1), k(2), k(3), and k(4),are vectors containing intermediate Runge-Kutta results– f is a vector containing the derivatives– k(1)= hfn= hf(xn, y(n))– k(2)= hf(xn+ h/2, y(n) + k(1)/2)– k(3)= hf(xn+ h/2, y(n) + k(2)/2)– k(4)= hf(xn+ h, y(n) + k(3))– y(n+1)= (k(1) + 2k(2) + 2k(3)+ k(4))/69ODE System by RK4• dy/dx = -y + z and dz/dx = y – z with y(0) = 1 and z(0) = -1 with h = .1• First compute k1for y and z•k(1)y= h[-y + z] = 0.1[-1 + (-1)] = -.2•k(1)z= h[y - z] = 0.1[1 - (-1)] = .2• Now we can compute k2for y and z•k(2)y= h[-(y+ k(1)y/2) + z + k(1)z/2] = 0.1[ -(1 + -0.2/2) + (-1 + .2/2)] = -.18•k(2)z= h[(y+ k(1)y/2) –(z + k(1)z/2)] = 0.1[(1 + -0.2)/2 - (-1 + .2/2)] = .1810ODE System by RK4 II• Continue as on previous chart•k(3)y= h[-(y+ k(2)y/2) + z + k(2)z/2] = 0.1[ -(1 + -0.18/2) + (-1 + .18/2)] = -.182•k(3)z= h[(y+ k(2)y/2) –(z + k(2)z/2)] = 0.1[(1 + -0.18)/2 - (-1 + .18/2)] = .182•k(4)y= h[-(y+ k(3)y) + z + k(3)z] = 0.1[ -(1 + -0.182) + (-1 + .182)] = -.1636•k(4)z= h[(y+ k(3)y) –(z + k(3)z)] = 0.1[ (1 + -0.182) - (-1 + .182)] = .163611ODE System by RK4 III•yi+1= yi+ (k(1)y+ 2k(2)y+ 2k(3)y+ k(4)y)/6 = 1 + [ (–.2) + 2(–.18) + 2(–.182) + (–.1636)]/6 = .8187•zi+1= zi+ (k(1)z+ 2k(2)z+ 2k(3)z+ k(4)z)/6 = –1 + [(.2) + 2(.18) + 2(.182) + (.1636)]/6 = –.8187• Continue in this fashion until desired final x value is reached• No x dependence for f in this example12Numerical Software for ODEs• Usually written to solve a system of N equations, but will work for N = 1• User has to code a subroutine or function to compute the f array– Input variables are x and y; f is output – Some codes have one dimensional parameter array to pass additional information from main program into the function that computes derivatives313Derivative Subroutine Example• Visual Basic code for system of ODEs at right is shown belowSub fsub( x as Double, y() as Double, _ f() as Double )f(1) = -y(1) + Sqr(y(2)) + y(3)*Exp(2*x)f(2) = -2 * y(1)^2f(3) = -3 * y(1) * y(2)End Sub2132122321132 yydxdyydxdyeyyydxdyx−=−=++−=14Derivative Subroutine Example• Fortran code for system of ODEsat right is shown belowsubroutine fsub( x, y, f )real(KIND=8) x, y(:), f(:)f(1) = -y(1) + sqrt(y(2)) +y(3)*exp(2*x)f(2) = -2 * y(1)**2f(3) = -3 * y(1) * y(2)end subroutine fsub2132122321132 yydxdyydxdyeyyydxdyx−=−=++−=15Derivative Subroutine Example• C++ code for system of ODEsat right is shown belowvoid fsub(double x, double y[], double f[]){f[1] = -y[1] + sqrt(y[2]) + y[3]*exp(2*x);f[2] = -2 * y[1] * y[1];f[3] = -3 * y[1] * y[2];}2132122321132 yydxdyydxdyeyyydxdyx−=−=++−=16Stiff Systems of Equations• Several definitions• Basic problem is that there are several length or time constants (eigenvalues) in the system– If one is negative and large in magnitude compared to others, this will set stability– Such terms quickly drop to zero and do not affect physical solution– However, they force small h for stability17Single Stiff Equation• Hoffman uses nonhomogenous equation with different length (time) scales • dy/dx = f(x,y) = -a[y – F(x)] + dF/dx• y = [y0– F(0)]e-ax+ F(x)– Solution details on next two charts• F(x) = cx + b: y = [y0–b] e-ax+ cx + b–y/y0= (1 – b/y0)e-ax+ cx/y0+ b/y0–y/y0= (1 – b/y0)e-ax+ ax(c/ay0) + b/y0–y/y0= g(ax); y/y0 → e-axas x → 0 and y/y0 → b/y0+ (c/ay0)(ax) as x →∞18Stiff Equation Solution• dy/dx = f(x,y) = -a[y – F(x)] + dF/dx solu-tion from formula for a first-order ODEdxdFaFxgaxfxgyxfdxdy+===+ )()()()(⎥⎦⎤⎢⎣⎡⎟⎠⎞⎜⎝⎛∫+∫=∫−dxexgCeydxxfdxxf )()()(⎥⎦⎤⎢⎣⎡⎟⎠⎞⎜⎝⎛∫+∫=∫−dxexgCeyadxadx)([]∫+=−dxexgCeyaxax)(419Stiff Equation Solution IIdxdxdFaFedxexgaxax⎟⎠⎞⎜⎝⎛+=∫∫)(∫∫∫=+= aFdxedFeaFdxeaxaxax[][])()( xFCeFeCedxexgCeyaxaxaxaxax+=+=+=−−−∫[])()0()0(00xFeFyysoFyCax+−=−=−By parts ∫udF = uF - ∫Fdu∫−+ )( dxaeFFeaxaxFeax=20Gear’s Method• Implicit algorithm• General algorithm shown below– k is global order of the method– Coefficients α-j, γ, and β in Hoffman Table 7.19, page 407, for various k values∑=−−+++=kjjnjnnyhfy011αγγβ• Usual problem of