First Order Linear Differential Equations16



3.8 Numerical Solutions to Systems of Differential Equations

A linear or nonlinear first order differential equation can always be solved numerically. Consider the following differential equation

[pic] = f(x, y) (3.8-1)

with initial condition x = x0, y = y0. The solution to the equation (3.8-1) is plotted in Figure 3.8-1.

[pic]

Figure 3.8-1. Numerical evaluation of [pic] = f(x, y)

The dependent variable y at x = x0 + h = x1 can be approximated by y = y0 + (y where

(y =h f(x0, y0)

The numerical solution y1 will approach the analytical solution as h approaches zero. The process can be repeated until the final value of the independent variable x is obtained. This procedure is the well-known Euler’s method to evaluate a first order differential equation numerically. The value at the end of each interval or step size h is evaluated from the values and slope at the beginning of the interval. This method only requires one slope evaluation therefore it is a first order method. For a given step size, the accuracy of the solution can be improved by evaluating the slope f(x, y) more than one. Each successive slope evaluation is based on the improved value of the dependent variable from the previous slope. The slope can be evaluated at any point between the interval h. However, it is normally evaluated at the beginning, the midpoint, and/or the end of the interval.

Runge-Kutta method

The Runge-Kutta Method achieves the accuracy of a Taylor series approach without requiring the calculation of higher derivatives.

Fourth order Runge-Kutta method for [pic] = f(x,y)

yn+1 = yn + [pic](k1 + 2k2 + 2k3 + k4)

where

k1 = h(f(xn, yn) k2 = h(f(xn + 0.5h, yn + 0.5k1)

k3 = h(f(xn + 0.5h, yn + 0.5k2) k4 = h(f(xn + h, yn + k3)

k2/h and k3/h are the slopes at the mid point. Each of the k/h represents a slope over the given interval. k1/h is the slope at the beginning of the interval and k4/h is the slope at the end of interval. Therefore the sum (k1 + 2k2 + 2k3 + k4) / (6h) can be interpreted as a weighted average of the slope. Due to the recurrence of the k parameters, that is k1 appears in the equation for k2, which appears in the equation for k3, and so forth, the Runge-Kutta method is very efficient for computer calculations.

EX: Solve [pic] = x + y ; y(0) = 1, h = 0.1

k1 = 0.1*(0 + 1) = 0.1*1 => y = y(0) + 0.5k1 = 1 + 0.5*0.1*1 = 1.05

k2 = 0.1*(0.05 + 1.05) = 0.1*1.10 => y = y(0) + 0.5k2 = 1 + 0.5*0.1*1.10 = 1.055

k3 = 0.1*(0.05 + 1.055) = 0.1*1.105 => y = y(0) + k3 = 1 + 1*1.105 = 1.1105

k4 = 0.1*(0.1 + 1.1105) = 0.1*1.2105

y(0.1) = 1 + 0.1*(1 + 2*1.1 + 2*1.105 + 1.2105)/6 = 1.11034

k1 = 0.1*(0.1 + 1.11034) = 0.1*1.21034

=> y = y(0.1) + 0.5k1 = 1.11034 + 0.5*0.1*1.21034 = 1.170857

k2 = 0.1*(0.15 + 1.170857) = 0.1*1.32857

=> y = y(0.1) + 0.5k2 = 1.11034 + 0.5*0.1*1.32857 = 1.17638285

Fourth-order Runge-Kutta method for system of first order differential equations

[pic]= f1(x, y1, y2, …, ym)

[pic]= f2(x, y1, y2, …, ym) ......

[pic]= fm(x, y1, y2, …, ym)

[pic] = [pic] + [pic] (k1,i + 2k2,i + 2k3,i + k4,i) , where i = 1, 2, ..., m and

k1,i = h*fi(xn, [pic], …, [pic]) k2,i = h*fi(xn + [pic], [pic] + [pic], …, [pic]+ [pic])

k3,i = h*fi(xn + [pic], [pic] + [pic], …, [pic]+ [pic])

k4,i = h*fi(xn + h, [pic] + k3,i, …, [pic]+ k3,m)

The idea behind the solution to a system of differential equations is similar to the solution of a single differential equation. All the k’s values divided by step size h are just the slopes of the curves evaluated at the appropriate x and yi.

EX: Solve [pic]= y1y2 + x , y1(0) = 1

[pic] = xy2 + y1 , y2(0) = -1

using fourth order Runge-Kutta method with step size h = 0.1

Solution

x = xn = 0, y1 = 1, y2 = -1

k1,1 = 0.1*(y1y2 + x) = 0.1*[ (1)(-1) + 0] = -0.1

k1,2 = 0.1*(xy2 + y1) = 0.1*[ (0)(-1) + 1] = 0.1

x = xn + 0.5h = 0 + (0.5)(.1) = 0.05

y1 = y1(0) + 0.5k1,1 = 1 + 0.5(-0.1) = 0.95

y2 = y2(0) + 0.5k1,2 = -1 + 0.5( 0.1) = -0.95

k2,1 = 0.1*(y1y2 + x) = 0.1*[ (0.95)(-0.95) + 0.05] = -0.08525

k2,2 = 0.1*(xy2 + y1) = 0.1*[ (0.05)(-0.95) + 0.95] = 0.9025

x = xn + 0.5h = 0 + (0.5)(.1) = 0.05

y1 = y1(0) + 0.5k2,1 = 1 + 0.5(-0.0852) = 0.9574

y2 = y2(0) + 0.5k2,2 = -1 + 0.5( 0.0902) = -0.9549

k3,1 = 0.1*(y1y2 + x) = 0.1*[ (0.9574)(-0.9549) + 0.05] = -0.0864

k3,2 = 0.1*(xy2 + y1) = 0.1*[ (0.05)(-0.9549) + 0.9574] = 0.0909

x = xn + h = 0 + .1 = 0.1

y1 = y1(0) + k3,1 = 1 + (-0.0864) = 0.9136

y2 = y2(0) + k3,2 = -1 + ( 0.0909) = -0.9091

k4,1 = 0.1*(y1y2 + x) = 0.1*[ (0.9136)(-0.9091) + 0.1] = -0.0730

k4,2 = 0.1*(xy2 + y1) = 0.1*[ (0.1)(-0.9091) + 0.9136] = 0.0823

y1(0.1) = y1(0) + (k1,1 + 2k2,1 + 2k3,1 + k4,1)/6

y1(0.1) = 1 + [(-0.1) + 2(-0.08525) + 2(-0.0864) + (-0.0730)]/6 = 0.9139

y2(0.1) = y2(0) + (k1,2 + 2k2,2 + 2k3,2 + k4,2)/6

y2(0.1) = -1 + [(0.1) + 2(0.09025) + 2(0.0909) + (0.0823)]/6 = -0.9092

x = xn = 0.1, y1 = 0.9139, y2 = -0.9092

k1,1 = 0.1*(y1y2 + x) = 0.1*[ (0.9139)(-0.9092) + 0.1] = -0.07309

k1,2 = 0.1*(xy2 + y1) = 0.1*[ (0.1)(-0.9092) + 0.9139] = 0.082298

Example 3.8-1. Solve the following first order system for y1 and y2 at x = 1.

[pic]= y1y2 + x , y1(0) = 1

[pic] = xy2 + y1 , y2(0) = -1

using fourth order Runge-Kutta method with step size h = 0.1

Solution

The MATLAB routines ode23 and ode45 can be used to solve the system. A MATLAB function must be created to evaluate the slopes as a column vector. The function name in this example is exode(x, y) which must be saved first in the hard drive with the same name exode.m.

-------------------------------

% Example 3.8-1 function to evaluate the slopes

function y12 = exode(x,y)

y12(1,1)=y(1)*y(2)+x;

y12(2,1)=x*y(2)+y(1);

-------------------------------

The command ode23 or ode45 is then evaluated from the command windows. MATLAB will set the step size to achieve a preset accuracy that can be changed by user. We will use both ode23 and ode45.

[x,y]=ode23('exode',[0; 1],[1; -1])

|x = |y = |

|0 |1.0000 -1.0000 |

|0.0800 |0.9290 -0.9260 |

|0.1800 |0.8628 -0.8480 |

|0.2800 |0.8174 -0.7828 |

|0.3800 |0.7899 -0.7275 |

|0.4800 |0.7779 -0.6794 |

|0.5800 |0.7797 -0.6364 |

|0.6800 |0.7943 -0.5966 |

|0.7800 |0.8207 -0.5581 |

|0.8800 |0.8585 -0.5189 |

|0.9800 |0.9076 -0.4770 |

|1.0000 |0.9188 -0.4681 |

The independent variable can also be specified at certain locations between the initial and final values and MATLAB will provide the dependent value at these locations. However, the step size h is still controlled by the error tolerance.

>> xspan=0:.1:1;

>> [x,y]=ode23('exode',xspan,[1 -1])

|x = |y = |

|0 |1.0000 -1.0000 |

|0.1000 |0.9139 -0.9092 |

|0.2000 |0.8522 -0.8341 |

|0.3000 |0.8106 -0.7711 |

|0.4000 |0.7863 -0.7173 |

|0.5000 |0.7772 -0.6704 |

|0.6000 |0.7816 -0.6283 |

|0.7000 |0.7986 -0.5889 |

|0.8000 |0.8274 -0.5503 |

|0.9000 |0.8675 -0.5108 |

|1.0000 |0.9188 -0.4681 |

>> [x,y]=ode45('exode',xspan,[1 -1])

|x = |y = |

|0 |1.0000 -1.0000 |

|0.1000 |0.9139 -0.9092 |

|0.2000 |0.8522 -0.8341 |

|0.3000 |0.8106 -0.7711 |

|0.4000 |0.7863 -0.7174 |

|0.5000 |0.7772 -0.6705 |

|0.6000 |0.7817 -0.6283 |

|0.7000 |0.7987 -0.5889 |

|0.8000 |0.8274 -0.5504 |

|0.9000 |0.8675 -0.5108 |

|1.0000 |0.9188 -0.4681 |

................
................

In order to avoid copyright disputes, this page is only a partial summary.

Google Online Preview   Download