Differential Equations Numerical Solution of ODEs in Matlab -- Lab Exercises

Numerical Solution of ODEs in Matlab -- Lab Exercises

The Euler Method

Getting Started -- Do the following:

Download a sequence of m-files from the website for this Differential Equations course:

http://www.eng.uml.edu/Dept/Chemical/onlinec/onlinec.htm

Get the following files and be sure to store them in a known location (a:\ or c:\temp):

Open Matlab and change the default directory to the desired location (a:\ or c:\temp).

Run the Matlab script file ndemo1.m by simply typing the name of the file at the Matlab prompt. This program is interactive in that it asks the user to input the number of integration steps to use over a given domain, and it does this repeatedly until the user types in a value of zero for the number of steps. Try various values, N = 5, N = 10, and N = 100, and you will see that the numerical solution of the ODE gets closer and closer to the exact solution as N increases or as the step size,

, decreases.

Well that was pretty easy -- you have just solved the IVP defined by

using the Euler Method for discretizing the ODE. You have also seen that "as decreases, the numerical solution becomes closer to the exact solution." Let’s now dissect these programs a little just to make sure you know what has happened.

Now What Did You Just Do? -- Try the following tasks:

Part I -- The function file to evaluate y’ = f(x,y)

Open the function file eqn1.m in the Matlab editor. The input arguments to this file are specific values of x and y. Knowing x and y, we can evaluate the function f(x,y) which is the derivative of y(x) at the point [x,y]. Thus, the goal of this routine is to return y’ = f(x,y) for the input x,y pair. In this case we give the result, y’, a variable name called yp and pass this variable back to the calling program as part of the function’s output argument list.

Following this discussion, let’s create a function file to handle the first problem for your HW assignment. To do this we can simply edit the current file by making the changes as indicated in the following table:

	Example Problem	HW Problem #1
Differential equation
Matlab syntax in function	*yp = exp(-2x)-y;**	yp = -y;

Now save this new file as p1_eqns.m to refer to the Prob. #1 equations (be sure to hit the Save As option in the File menu in the editor and be sure you know where the file is being saved).

Part II -- The Euler Method (euler.m and euler1.m)

Open the function file euler.m in the editor. This file implements a general form of the Euler Method. However, it is more complicated than necessary at this time. The added complexity gives this routine greater flexibility for later use in this course (i.e. it will handle a system of first-order ODEs). At this point we don’t need all this detail, so a single-equation version, called euler1.m, has been generated for discussion here. Open euler1.m and briefly compare it to euler.m. You should see that the two files are similar, but that euler1.m is slightly less complicated.

Now let’s focus on the euler1.m file. A listing of euler1.m is shown below:

%
%   EULER1.M   Use Euler Method for Numerical Solution of IVPs
%
%   This function solves a single first-order ODE via Euler's Method
%
%   The inputs to the function are:
%      fxy = string variable with the name of the function file containing f(x,y)
%      xo,xf = initial and final values of the independent variable (scalars)
%      yo = initial value of dependent variable at xo 
%      N = number of intervals to use between xo and xf
%   The outputs to the function are:
%      X = vector containing values of the independent variable
%      Y = the estimated dependent variable at each value of the independent variable
%

      function [X,Y] = euler1(fxy,xo,xf,yo,N)
%
%   compute step size and size of output variables
      if N < 2   N = 2;   end  % set minimum number for N
      h = (xf-xo)/N;           % step size
      X = zeros(N+1,1);        % initialize independent variable
      Y = zeros(N+1,1);        % initialize dependent variables
%
%   set initial conditions
      X(1) = xo;   Y(1) = yo;
%
%   begin computational loop
      for i = 1:N
        k1 = h*feval(fxy,X(i),Y(i));   % evaluate function
        Y(i+1) = Y(i) + k1;            % increment dependent variable
        X(i+1) = X(i) + h;             % increment independent variable
      end
%
%   end of function
%

The first part of this file computes the step size, h, and does some initialization. The final few lines of code contain the basic Euler algorithm:

Loop over number of time steps	Loop over all i
Evaluate function
Increment dependent variable
Increment independent variable
End loop	End loop

Note that Matlab’s feval function is used to implement the actual function evaluation. This is needed since we are using the name of a generic function, fxy, in the current routine. The command feval(fxy,x,y) is essentially equivalent to eqn1(x,y) or p1_eqns(x,y), where the local variable, fxy, has been equivalenced to the name of the ODE function file, eqn1 or p1_eqns, via the input argument list in these examples.

Now let’s test drive euler1.m. For example, let’s solve the base problem from above with the following commands:

xo = 0; xf = 0.5; yo = 1; N = 10;
[x,y] = euler1('eqn1',xo,xf,yo,N); plot(x,y)

Pretty easy, don’t you think? This should look similar to the plots generated above (except we have not included the exact solution). Using the same basic scheme, you can solve the differential equation associated with your HW assignment. For example, Prob. #1 becomes

and we already have the ODE function built into function file p1_eqns.m. Thus, solution of this problem can be accomplished by

xo = 0; xf = 0.5; yo = 2; N = 10;
[x,y] = euler1('p1_eqns',xo,xf,yo,N); plot(x,y)

This is pretty crude but it contains the basics. Note that euler.m could be used if desired.

Part III -- Putting it all together into a Matlab program -- ndemo1.m

The program ndemo1.m simply elaborates on the above basic solution scheme. It includes the evaluation of the exact solution for comparison purposes and for studying the sensitivity of the solution to the choice of step size. It also includes some intermediate results that are tabulated to help compare with the hand calculations you were asked to do. Finally, it enhances the plots a little with a title, x-axis label, etc.. Putting everything into a program let’s you do more things -- these all could be done interactively but it could become pretty tedious. The program also can be used as formal documentation of exactly what computations were performed. Following this guide, you should modify ndemo1.m, as needed, to solve the requested HW problems.

You are now an expert with using the Euler Method for solving simple IVPs.
Congratulations!!!

92.236 Matlab Lab Exercises by Dr. J. R. White, UMass-Lowell (Sept. 2000).

Return to Online Courses