A modification to the regula falsi method uses two starting estimates that must bracket the root. The algorithm is as follows:

**1.a.**Given two starting estimates, x1 and x2 that bracket the root, compute the mid point, xm = (x1 + x2)/2.

**b.**Compute the corresponding function evaluations for these three points: y1 = f(x1), y2 = f(x2) and ym = f(xm).

**c.** Compute s = If s = 0, then set x = xm and stop. Otherwise continue.

**d.** Update the estimate

**e.** Tighten the bracket by keeping xn and one of x1, x2 or xm with:

If xm and xn bracket the root, set x1 = xm and x2 = xn

Else, if x1 and xn bracket the root, set x2 = xn

Else, it must be the case that xn and x2 bracket the root, so set x1 = xn.

**f.** Continue by returning to step 1.

Write a Matlab root finding routine using the following algorithm. Test your algorithm on a variety of nonlinear functions. How does it compare to other algorithms?

ODEs

**2.** A simple model of the spread of disease gives

where P(t) represents the number of individuals in the population who are infected, and C is the constant size of the total population.

**(a)** Show that the solution of this differential equation is the logistic function.

**(b)** Suppose that the parameter k fluctuates (perhaps because the population is more susceptible during certain seasons). Solve the modified problem

with P(0) = 10, k = 2, C = 200.

**3.** Solve the Bessel equation of order one,

starting with y(1) = 1 and y 0 (1) = 0 on 1 ≤ x ≤ 4

**4.** The motion of one body around another, such as a comet orbiting around the sun, can be described by the system

where the distance r =

With distance in AU (1 AU = 1.496 × 1011 m) and time in years, we have K ≈ 40 for an object rotating around the sun.

Take as initial conditions x(0) = 1, x0 (0) = 0 and y(0) = 0, y0 (0) = 2, and solve for 0 ≤ t ≤ 4.

Investigate the effects of different values of y 0 (0).

**5.** The equations for the deflection y and rotation z or a simply supported beam with a uniformly distributed load of intensity 9 · 103 N/ft and bending moment M(x) = 10x−x² can be expressed as

where E is the modulus of elasticity, and I is the moment of inertia of the cross section of the beam.

Taking EI = 1.6 · 107 N and y(0) = 0 and z(0) = −0.02, find y and z for 0 < x < 10.

**6.** Numerically compare the following 4th order Runga-Kutta-Gill method with the classical Runge-Kutta method.

## Reviews

There are no reviews yet.