the Consider simple differential equation The solution is y e
10t
y ' 10 y
y (0) 1
Consider now the solution of this equation with Euler's method.
y (t t ) y(t ) y '(t )t y (t ) 10 y (t )t y (t )(1 10t )
If we take a large t the solution will be inaccurate, and if t is larger than 0.2 it will be unstable. If we are worried about stability, we can use the implicit Euler method
y (t t ) y (t ) y '(t t )t y (t ) 10 y (t t )t y (t t )(1 10t ) y (t ) y (t t ) y (t ) / (1 10t )
In terms of accuracy the two methods are comparable. Consider for example taking t 0.1 With Euler's method (explicit) we obtain
y (0.1) 1 10t 0
With Euler's implicit method we get
y (0.1) 1/ (1 10t ) 0.5
The exact solution is e
1
0.3679
Heun's method takes the average of the derivative at the beginning and at the predicted point. That is, after we predict that y(0.1)=0, the derivative there is also zero. So we recalculate with a derivative that is the average of 10 (the value at t=0) and zero, so
yHeun (0.1) 1 5 0.1 0.5
This is an example of a predictorcorrector method. We make a prediction with Euler, and we correct it with Heun's. We can take one additional correction noting that at 0.5 the derivative is 5, so the new average is 7.5, which will predict y(0.1)=0.25. We can also solve it with the ode45 function >> yfun=@(t,y) 10*y
yfun = @(t,y) 10*y >> [tout,yout]=ode45(yfun,[0 1],1); >> tout' ans = Columns 1 through 10 0 0.0050 0.0100 0.0151 0.0201 0.0406 0.0612 0.0817 0.1023 0.1219 Columns through 11 20 0.1415 0.1611 0.1807 0.2004 0.2200 0.2397 0.2593 0.2790 0.2987 0.3183 Columns 21 through 30 0.3380 0.3576 0.3773 0.3970 0.4166 0.4363 0.4559 0.4756 0.4953 0.5149 Columns 31 through 40 0.5346 0.5542 0.5739 0.5936 0.6132 0.6329 0.6525 0.6722 0.6918 0.7115 Columns 41 through 50 0.7312 0.7508 0.7705 0.7901 0.8098 0.8311 0.8524 0.8737 0.8950 0.9213 Columns 51 through 53 0.9475 0.9738 1.0000 >> format short e >> yout' ans = Columns 1 through 7 1.0000e+000 9.5100e001 9.0441e001 8.6009e001 8.1795e001 6.6588e001 5.4211e001 Columns 8 through 14
4.4161e001 3.5979e001 2.9565e001 2.4296e001 1.9975e001 1.6424e001 1.3489e001 Columns 15 through 21 1.1079e001 9.1041e002 7.4818e002 6.1451e002 5.0474e002 4.1476e002 3.4086e002 Columns 22 through 28 2.7996e002 2.2995e002 1.8896e002 1.5530e002 1.2755e002 1.0477e002 8.6090e003 Columns 29 through 35 7.0751e003 5.8110e003 4.7730e003 3.9222e003 3.2234e003 2.6475e003 2.1745e003 Columns 36 through 42 1.7869e003 1.4685e003 1.2062e003 9.9070e004 8.1410e004 6.6905e004 5.4952e004 Columns 43 through 49 4.5135e004 3.7090e004 3.0481e004 2.4624e004 1.9893e004 1.6083e004 1.3005e004 Columns 50 through 53 9.9907e005 7.6759e005 5.9137e005 4.5600e005 We can try to shorten the calculation with lower accuracy requirements >> options=odeset('RelTol',1,'AbsTol',1,'InitialStep',0.1); However, this does not appear to help.