Question

[7] 1. Consider the initial value problem (IVP)

y′(t) = −y(t), y(0) = 1

The solution to this IVP is

y(t) = e−t

[1] i) Implement Euler’s method and generate an approximate solution of this IVP over the interval [0,2], using stepsize h = 0.1. (The Google sheet posted on LEARN is set up to carry out precisely this task.) Report the resulting approximation of the value y(2).

[1] ii) Repeat part (ii), but use stepsize h = 0.05. Describe the process you followed to set up the calculation, and report the resulting approximation of the value y(2).

[2] iii) Determine the absolute error in the approximations generated in parts (i) and (ii). The stepsize is reduced by a factor of one half from (i) to (ii). By what factor is the error reduced (roughly)?

[1] iv) Using your result from part (iii), what would you conjecture the error to be if the simulation is repeated with stepsize h = 0.025?

[2] v) Attempt to generate an approximate solution of this IVP using Euler’s method over the interval [0, 10], using stepsize h = 2. Explain your results.

[4] 2. Consider the initial value problem (IVP)

y′(t) = 4e−ty(t) − [y(t)]2 y(0) = 2

This differential equation is neither separable nor linear, so we have no technique to solve it analytically (i.e. with paper and pencil). Implement Euler’s method and generate an approximate solution of this IVP over the interval [0, 3], using stepsize h = 0.1. Provide a description of your implementation of the algorithm and a graph showing your simulation.

[9] 3. Recall that Euler’s method for the differential equation y′(t) = f(t,y(t))

1

is derived from the approximation

by setting y′(t) = f(t,y(t)): and then solving for y(t + h).

y′(t) ≈ y(t + h) − y(t) h

f(t,y(t))≈ y(t+h)−y(t) h

Applying the same procedure, but starting with the approximation y′(t) ≈ y(t) − y(t − h)

h

we arrive at an update rule of the form

y(t) = y(t − h) + hf (t, y(t)) This update rule, which is more typically written as

y(t + h) = y(t) + hf (t + h, y(t + h))

is called the implicit Euler method. This update rule does not appear to be useful: it does not provide a formula for the ‘next’ value y(t + h). Instead, it provides only an equation, which must be solved to determine the value of y(t + h). As it turns out, the implict Euler method has several advantages over the regular (explicit) Euler method, and so this extra effort is often worth the trouble. In general, this equation must be solved numerically, by, e.g. applying Newton’s method (at each step of the simulation).

[3] i) Derive the implicit Euler update rule for the differential equation: y′(t) = 4e−ty(t) − [y(t)]2 y(0) = 2

In this case Newton’s method is not needed, because the update equation can be solved explicitly for y(t+h). You will need to solve a quadratic equation to determine y(t + h) as a function of y(t). You’ll need to use the fact that the solution of this IVP only takes positive values.

[3] ii) Implement you update law from part (ii) and generate an approximate solution of this IVP over the interval [0, 3], using stepsize h = 0.1. Provide a description of your implementation of the algorithm and a graph showing your simulation. Comment on whether your results agree with the simulation in Question 2.

[3] iii) Simulations generated by either the original (explicit) Euler method or the implicit Euler method accurately estimate solution behaviour for very small values of stepsize h. However, a significant advantage of the implicit approach is that it can often provide reasonably accurate simulations even when the stepsize is chosen fairly large, which can be used for efficient computation. As an example, consider the solution of the IVP in Question 2 (and part (ii) of this question). The solution values rise, then fall and converge toward zero. This trend is clear when examining the solution over the time interval [0, 20]. Compare the efficiency of simulations that accurately capture the long-term solution trend by simulating the solution with both methods. That is, compare the stepsize choices that provide reasonable simulations for the two methods. Report on your findings.

Answer 1

%%% (a)

for h=0.1;

t y(t)
0.000000 1.000000
0.100000 0.900000
0.200000 0.810000
0.300000 0.729000
0.400000 0.656100
0.500000 0.590490
0.600000 0.531441
0.700000 0.478297
0.800000 0.430467
0.900000 0.387420
1.000000 0.348678
1.100000 0.313811
1.200000 0.282430
1.300000 0.254187
1.400000 0.228768
1.500000 0.205891
1.600000 0.185302
1.700000 0.166772
1.800000 0.150095
1.900000 0.135085
2.000000 0.121577

(ii) for h=0.05;

t y(t)
0.000000 1.000000
0.050000 0.950000
0.100000 0.902500
0.150000 0.857375
0.200000 0.814506
0.250000 0.773781
0.300000 0.735092
0.350000 0.698337
0.400000 0.663420
0.450000 0.630249
0.500000 0.598737
0.550000 0.568800
0.600000 0.540360
0.650000 0.513342
0.700000 0.487675
0.750000 0.463291
0.800000 0.440127
0.850000 0.418120
0.900000 0.397214
0.950000 0.377354
1.000000 0.358486
1.050000 0.340562
1.100000 0.323534
1.150000 0.307357
1.200000 0.291989
1.250000 0.277390
1.300000 0.263520
1.350000 0.250344
1.400000 0.237827
1.450000 0.225936
1.500000 0.214639
1.550000 0.203907
1.600000 0.193711
1.650000 0.184026
1.700000 0.174825
1.750000 0.166083
1.800000 0.157779
1.850000 0.149890
1.900000 0.142396
1.950000 0.135276
2.000000 0.128512

(iii)

y(2) actual value is exp(-2)=0.1353

error in (i) = abs(0.1353- 0.121577 ) = 0.0137

error in (ii) = abs(0.1353- 0.128512 ) = 0.0068

factor by which error reduced = 0.0137/0.0068=2.0147= 2

(iv)

thus for h=0.025 error will reduced by 2 = 0.0068/2= 0.0034