Leveraging the Power of Java and Matlab to Solve ODE's

Ordinary Differential Equations (ODE) are used to model a wide range of physical processes. An ODE is an equation containing a function of one independent variable and its ordinary derivatives. This paper presents the development and application of a practical teaching module introducing java programming techniques to electronics, computer, and bioengineering students before they encounter digital signal processing and its applications in junior and senior level courses. This paper will focus primarily on how to solve ODEs using Java and Matlab programming tools. There are two basic types of boundary condition categories for ODEs – initial value problems and two-point boundary value problems. Initial value problems are simpler to solve because you only have to integrate the ODE one time. The solution of a twopoint boundary value problem usually involves iterating between the values at the beginning and end of the range of integration. Runge-Kutta schemes are among the most commonly used techniques to solve initial-value problem ODEs. Matlab also presents several tools for modeling linear systems. These tools can be used to solve differential equations arising in such models, and to visualize the input-output relations. This paper attempts to describe how to use Java programming tool to solve initial value problems of ordinary differential equations (ODEs) using the Runge-Kutta scheme. It will also discuss how to represent initial value problems and demonstrate how to apply Matlab’s ODE solvers to such problems. It will also explain how to select a solver and how to specify solver options for efficient, customized execution. This paper provides an introduction to the Ordinary Differential Equations(ODEs). After a quick overview of selected numerical methods for solving differential equations using Matlab, we will briefly give an account of Euler and modified Euler methods for solving first order differential equations. This will be followed by numerical method for systems specially RungeKutta schemes and applications of second order differential equations in mechanical vibrations and electric circuits by leveraging the power of Java and Matlab. This paper will explain how this learning and teaching module is instrumental for progressive learning of students; the paper will also demonstrate how the numerical and integral algorithms are derived and computed through leverage of the java data structures. As a result, there will be a discussion concerning the comparison of Java and Matlab programming as well as students’ feedback. The result of this new approach is expected to strengthen the capacity and quality of our undergraduate degree programs and enhance overall student learning and satisfaction. Introduction to ODEs Differential equations are used to model a wide range of physical processes; technology students will use them in chemistry, biophysics, mechanics, thermodynamics, electronics, and P ge 24866.2 almost every other scientific and engineering discipline. An ODE is used to express the rate of change of one quantity with respect to another. One defining characteristic of an ODE is that its derivatives are a function of one independent variable. The order of a differential equation is defined as the order of the highest derivative appearing in the equation and ODE can be of any order. A general form of a first-order ODE can be written in the form dy/dx + p(x)y + q(x) + r = 0 where p(x) and q(x) are functions of x. This equation can be rewritten as shown below d/dx(y) +y p(x) = q(x) r where r is zero. A classical integrating factor method can be used for solving this linear differential equation of first order. The integrating factor is (exp)^∫p dx. Euler Method Graphical methods produce plots of solutions to first order differential equations of the form y’ = f(x,y), where the derivative appears on the left side of the equation. If an initial condition of the form y(x0) = y0 is also specified, then the only solution curve of of interest is y’ = f(x,y) the one that passes through the intial point (x0,y0). For the first-order initial-value problem the popular graphical method also known as Euler method can be used that satisfies the formula given below yn+1 = yn + hf(xn ,yn ) which can also be written as yn+1 = yn + h(y’n ), where the approximate solution at xn is designated by y(xn), or simply yn. The true solution at xn will be denoted by either Y(xn) or Yn. Note that once yn is known, equation y’ = f(x,y) can be used to obtain yn’ as yn’ = f(xn ,yn ) (1.0) Modified Euler’s Method: This is a simple predictor-corrector method that uses Euler’s method as the predictor and then uses the average value of y’ at both the left and right end points of the interval [xn, xn+1] (n= 0,1,2, ...) as the slope of the line element approximation to the solution over that interval. The resuting equations are: predictor: yn+1 = yn + h(y’n ) corrector: yn+1 = yn + h/2 For notational convenience, we designate the predicted value of yn+1 by pyn+1. Since y’n = f(xn ,yn ), it then follows that py'n+1 = f(xn+1 ,yn+1 ) (1.1) and the modified Euler method becomes predictor: pyn+1 = yn + h(y’n ) corrector: yn+1 = yn + h/2(py’n+1 + y’n) (1.2) P ge 24866.3 Example 1 Use the modified Euler’s method to solve y’ – y + x = 0; y(0) =2 on the interval[0,1] with h= 0.1 Using Matlab dsolve function we get >> dsolve('Dy = y x', 'y(0)= 2','x') the solution y = x + exp(x) + 1 Euler’s modified Numerical Method Here f(x,y) = y – x, and y0 =2. From equation (1.0), we have y0 = f(0,2) = 2 – 0 = 2 Using equations (1.1) and (1.2), we compute For n = 0, x1 = 0.1 py1 = y0 + h y0’ = 2 + 0.1(2) = 2.2 py1’ = f (x1,py1) = f(0.1,2.2) = 2.2 – 0.1 = 2.1 y1 = y0 + h/2(py1’ + y0’) = 2 + 0.05(2.1 + 2) = 2.205 y1’ = f (x1,y1) = f (0.1, 2.205) = 2.205 0.1 = 2.105 It can be shown that For n = 1, x2 = 0.2 py2 = y1 + h y1’ = 2.4155 py2’ = f (x2,py2) = f(0.2,2.4155) = 2.2155 y2 = y1 + h/2(py2’ + y1’) = 2.421025 y2’ = f (x2,y2) = f (0.2, 2.421025) = 2.221025 and so on Instead of computing xn , yn, and true solution etc. at different points, the following Matlab script can be used to obtain the Modified Euler solution as depicted in Figure below. % nsteps = 150 Final_Time=10.0; stepsize=Final_Time / nsteps; clear x y exactSolution y(1)=2.205; x(1)=0.1; exactSolution(1)=4.718281 8; for k=1:nsteps x(k+1)=x(k) + stepsize; y(k+1)=y(k) + stepsize*((y(k) x(k))); exactSolution(k+1) = exp(x(k+1)) + x(k+1) + 1; end plot(x,y); % default line color is blue title(‘Modified Euler Method solution of yprime = y – x, y(0) = 2’) hold on plot(x,exactSolution,'g'); % g for green line legend('Euler solution','Exact P ge 24866.4 solution') hold off error=norm(y-exactSolution)/norm(exactSolution); A general form of a second-order ODE is shown as follows: d 2 y/dx 2 + p(x)dy/dx + q(x)y + r(x) + s = 0 (1.2) Any high order ODE can be expressed as a coupled set of first -order differential equations. For example the second-order ODE given in equation (1.2) can be reduced to a coupled set of two first-order differential equations. d/dx(dy/dx) = p(x)dy/dx – q(x)y – r(x) – s (1.3) d/dx(y) = dy/dx Java’s ODE Class We will use and demonstrate a class named ODESolver that will define a number of methods 2 used to solve ODEs and also subclasses that can be used to represent modeling real world applications. The ODE class will declare a field to store the number of first order equations. The other declared field will store the number of free variables (those that are not specified by boundary P ge 24866.5 conditions at the beginning of the integration range). For initial value problems, the number of free variables will be zero. Two-point boundary problems will have one or more free variables. The ODE class will declare methods to set the conditions at the start of the integration range and to compute the error at the end. These methods are ODE-specific, but since the ODE class represents a generic ODE, they will be implemented as stubs, and the required ODE subclasses will override those methods according to their needs. The Runge-Kutta family of methods are step-wise integration algorithms that are expected to give good results as long as very high accuracy is not required. Starting from an initial condition, the ODE is solved at discrete steps over the desired integration range. Consider a simple first order differential equation which can be integrated in a step-wise manner to solve for the dependent variable y. dy/dx = f(x,y) (1.4) dy = dx f(x,y) Replacing the derivative by its delta form, we can write Δy = yn+1 yn = Δx f(x,y) (1.5) Equation (1.5) is the general form of the equation that is solved. Using Taylor series expansion on the Euler’s method, we find that it is first order accurate in Δx and is only useful for linear or nearly linear functions. If the slope of the f(x,y) curve changes between changes xn and xn+1, the Euler’s method will compute an incorrect value for yn+1. This is why we choose Runge-Kutta method since it performs a successive approximation of yn+1 by evaluating the f(x,y) function at different locations between xn and xn+1. There are numerous Runge-Kutta schemes of various orders of accuracy, we will use one of the popular fourth-order algorithm accurate in Δx. This algorithm consists of five steps, four successive approximations of Δy and a fifth step that computes yn+1 based on a linear combination of the successive approximations. The fourth order Runge-Kutta method solution process steps are as follows: 1. Compute Δy1 using Euler’s method Δy1 = Δx f (xn, yn) 2. Compute Δy2 by evaluating f(x,y) at (xn + 1⁄2 Δx, yn + 1⁄2 Δy1) Δy2 = Δx f (xn + 1⁄2 Δx, yn + 1⁄2 Δy1) 3. Compute Δy3 by evaluating f(x,y) at (xn + 1⁄2 Δx, yn + 1⁄2 Δy2) Δy3 = Δx f (xn + 1⁄2 Δx, yn + 1⁄2 Δy2) 4. Compute Δy4 by evaluating f(x,y) at f (xn + Δx, yn + Δy3) P ge 24866.6 Δy4 = Δx f (xn + Δx, yn + Δy3) 5. Compute yn+1 using a linear combination of Δy1 thru Δy4 such as Δy1/6 + Δy2/3 + Δy3/3 + Δy