MATLAB With ODE-Modeling and Simulation-Lecture Handouts, Lecture notes of Mathematical Modeling and Simulation

Download MATLAB With ODE-Modeling and Simulation-Lecture Handouts and more Lecture notes Mathematical Modeling and Simulation in PDF only on Docsity!

unknown function is called an ordinary differential equation, abbreviated as ODE. The order of the equation is determined by the order of the highest derivative. For example, if the first derivative is the only derivative, the equation is called a first-order ODE. In the same way, if the highest derivative is second order, the equation is called a second-order ODE.

The problems of solving an ODE are classified into initial-value problems (IVP) and boundary- value problems (BVP), depending on how the conditions at the endpoints of the domain are spec- ified. All the conditions of an initial-value problem are specified at the initial point. On the other hand, the problem becomes a boundary-value problem if the conditions are needed for both initial and final points. The ODE in the time domain are initial-value problems, so all the conditions are specified at the initial time, such as t = 0 or x = 0. For notations, we use t or x as an independent variable. Some literatures use t as time for independent variable.

It is important to note that our focus here is on the practical use of numerical methods in order to solve some typical problems, not to present any consistent theoretical background. There are many excellent and exhaustive texts on these subjects that may be consulted. For example, we would recommend Edwards and Penny (2000) [2], Boyce and DiPrima (2001) [3], Coombes et al. (2000) [4], Van Loan (1997) [5], Nakamura (2002) [6], Moler (2004) [7], and Gilat (2004) [8].

2 Numerical methods

Numerical methods are commonly used for solving mathematical problems that are formulated in science and engineering where it is difficult or even impossible to obtain exact solutions. Only a limited number of differential equations can be solved analytically. Numerical methods, on the other hand, can give an approximate solution to (almost) any equation. An ordinary differential equation (ODE) is an equation that contains an independent variable, a dependent variable, and derivatives of the dependent variable. Literal implementation of this procedure results in Euler’s method, which is, however, not recommended for any practical use. There are other methods more sophisticated than Euler’s. Among them, there are three major types of practical numerical methods for solving initial value problems for ODEs: (i) Runge-Kutta methods, (ii) Burlirsch-Stoer method, and (iii) predictor-corrector methods. We will present these three approaches on another occasion. Now, we are interested to talk about Euler’s methods.

2.1 EULER methods

The Euler methods are simple methods of solving first-order ODE, particularly suitable for quick programming because of their great simplicity, although their accuracy is not high. Euler methods include three versions, namely,

• forward Euler method
• modified Euler method
• backward Euler method

2.2.1 Basic steps

The typical steps of Euler’s method are given below.

Step 1. define f (x, y)

Step 2. input initial values x 0 and y 0

Step 3. input step sizes h and number of steps n

Step 4. calculate x and y:

for i=1:n x=x+h y=y+hf(x,y) end

Step 5. output x and y

Step 6. end

2.2.2 Example

As an application, consider the following initial value problem

dy dx

x y , y(0) = 1 (6)

which was chosen because we know the analytical solution and we can use it for check. Its exact or analytical solution is found to be

y(x) =

√ x^2 + 1 (7)

Therefore, we will be able to compare the approximate solutions and the exact solution. Here we wish to approximate y(0.3) using the Euler’s methods with step sizes h = 0.1 and h = 0.05. We find by hand-calculation,

x 0 = 0 , x 1 = 0. 1 , x 2 = 0. 2 , x 3 = 0. 3 y 0 = 1 y 1 = y 0 + hf (x 0 , y 0 ) = y 0 + hx 0 /y 0 = 1 y 2 = y 1 + hx 1 /y 1 = 1. 01 y 3 = y 2 + hx 2 /y 2 = 1. 0298.

Since y(0.3) =

√ (0.3)^2 + 1 = 1.044030, we find that

Error =

|y(0.3) − y 3 | y(0.3)

× 100 = 1.36%

Similarly, for the step size h = 0.05, we find that the error is

Error =

|y(0.3) − y 6 | y(0.3)

× 100 = 0.67%

2.2.3 MATLAB codes

• Step 1: Create user-defined function files: euler forward.m, euler modified.m, and euler backward.m.

function [x,y]=euler_forward(f,xinit,yinit,xfinal,n) % Euler approximation for ODE initial value problem % Euler forward method % File prepared by David Houcque - Northwestern U. 5/11/

% Calculation of h from xinit, xfinal, and n h=(xfinal-xinit)/n;

% Initialization of x and y as column vectors x=[xinit zeros(1,n)]; y=[yinit zeros(1,n)];

% Calculation of x and y for i=1:n x(i+1)=x(i)+h; y(i+1)=y(i)+h*f(x(i),y(i)); end

end

function [x,y]=euler_modified(f,xinit,yinit,xfinal,n) % Euler approximation for ODE initial value problem % Euler modified method % File prepared by David Houcque - Northwestern U. - 5/11/

% Calculation of h from xinit, xfinal, and n h=(xfinal-xinit)/n;

% Initialization of x and y as column vectors x=[xinit zeros(1,n)]; y=[yinit zeros(1,n)];

% Call functions [x1,y1]=euler_forward(f,0,1,0.3,6); [x2,y2]=euler_modified(f,0,1,0.3,6); [x3,y3]=euler_backward(f,0,1,0.3,6);

% Plot plot(xe,ye,’k-’,x1,y1,’k-.’,x2,y2,’k:’,x3,y3,’k--’) xlabel(’x’) ylabel(’y’) legend(’Analytical’,’Forward’,’Modified’,’Backward’) axis([0 0.3 1 1.07])

% Estimate errors error1=[’Forward error: ’ num2str(-100(ye(end)-y1(end))/ye(end)) ’%’]; error2=[’Modified error: ’ num2str(-100(ye(end)-y2(end))/ye(end)) ’%’]; error3=[’Backward error: ’ num2str(-100*(ye(end)-y3(end))/ye(end)) ’%’];

error={error1;error2;error3}; text(0.04,1.06,error)

• Step 3: Compare the results. The calculated results are displayed in the graphical form below. Reasonably good results are obtained even for a moderately large step size and the approximation can be improved by decreasing the step size. According to the results (Figure 1) and Table 1, forward and backward approaches give identically the same results (less than 1% of error), while modified method give very good result when compared with the exact solution.

h Forward Modified Backward

Table 1: Comparison of exact solution with Euler methods

2.3 Using built-in function

MATLAB has several different functions (built-ins) for the numerical solution of ordinary differ- ential equations (ODE). In this section, however, we will present one of them. We will also give an example how to use it, instead of writing our own MATLAB codes as we did in the first part. The basic steps, previously defined, are still typically the same. These solvers can be used with the following syntax:

[x,y] = solver(@odefun,tspan,y0)

0 0.05 0.1 0.15 0.2 0.25 0.

1

x

y

Forward error: −0.6751% Modified error: 0.0003983% Backward error: 0.67244%

Analytical Forward Modified Backward

Figure 1: Comparison of exact solution with Euler methods

Solver Accuracy Description

ode45 Medium This should be the first solver you try

ode23 Low Less accurate than ode

ode113 Low to high For computationally intensive problems

ode15s Low to medium Use if ode45 failed

Table 2: Some MATLAB ODE solvers

ode45(’xdot’,tspan,y0)

Note the use of @xdot and ’xdot’. Use function handles to pass any function that defines quantities the MATLAB solver will compute, in particular for simple functions.

On the other hand, it is also important to remember that complicated differential equations should be written an M-file instead of using inline command or function handle.

Following is an example of an ordinary differential equation using a MATLAB ODE solver. First, let’s create a script file, called main2.m, as follows:

% Script file: main2.m % The RHS of the differential equation is defined as % a function handle. % File prepared by David Houcque - Northwestern U. - 5/11/

f=@(x,y) x./y;

% Calculate exact solution g=@(x) sqrt(x.^2+1); xe=[0:0.01:0.3]; ye=g(xe);

% Call function [x4,y4]=ode45(f,[0,0.3],1);

% Plot plot(xe,ye,’k-’,x4,y4,’k:’) xlabel(’x’) ylabel(’y’) legend(’Analytical’,’ode45’) axis([0 0.3 1 1.05])

Here, we use the same data as defined in the first part for Euler’s methods. The initial conditions and the time steps are the same as before.

The integration proceeds by steps, taken to the values specified in tspan. Note that the step size (the distance between consecutive elements of tspan) does not have to be uniform.

A plot comparing the computed y versus x is shown in Figure 2. According to the plot, the

0 0.05 0.1 0.15 0.2 0.25 0.

1

x

y

Analytical ode

Figure 2: Comparison of exact solution with ode45 solution

calculated results using the built-in function (ode45) give a very good result when compared with the analytical solution.

Additional on all of these solvers can be found in the online help. For addition information on numerical methods, we refer to Shampine (1994) [10] and Forsythe et al. (1977) [11].

3 Conclusion

The above plots show the results obtained from different algorithms. Consequently, we can see better that Runge-Kutta algorithm is even more accurate at large step size (h = 0.1) than Euler algorithms at small step size (h = 0.05). One can easily adapt these MATLAB codes as needed for a different type of problem. In using numerical procedure, such as Euler’s method, one must always keep in mind the question of whether the results are accurate enough to be useful. In the preceding examples, the accuracy of the numerical results could be ascertained directly by a comparison with the solution obtained analytically. Of course, usually the analytical solution is not always available to compare.