1
00:00:00,160 --> 00:00:04,980
In the current chapter 12 of the numerical methods course, we study efficient numerical

2
00:00:05,147 --> 00:00:09,847
methods for the integration of stiff initial value problems.

3
00:00:10,014 --> 00:00:16,694
In the previous section, we learned about the family of General Runge Kutta single step

4
00:00:16,861 --> 00:00:17,961
method.

5
00:00:18,128 --> 00:00:24,668
We found that this family comprises quite a few members who do an excellent job about

6
00:00:24,835 --> 00:00:27,755
efficiently dealing with stiff initial value problems.

7
00:00:27,922 --> 00:00:32,842
This means they are not subject to stability induced time step constraints.

8
00:00:33,009 --> 00:00:40,689
Rather they are unconditionally stable and the choice of the time step size is completely

9
00:00:40,856 --> 00:00:43,896
dictated by the accuracy requirements.

10
00:00:44,063 --> 00:00:45,223
However there is a catch.

11
00:00:45,390 --> 00:00:52,150
We have seen by model problem analysis that necessarily all these unconditionally stable

12
00:00:52,317 --> 00:00:55,877
Runge Kutta single step methods have to be implicit.

13
00:00:56,044 --> 00:01:02,224
This means that in every time step we have to solve the increment equations and those

14
00:01:02,391 --> 00:01:09,771
amount to well an ns times ns nonlinear system of equations and its solution can be expensive.

15
00:01:09,938 --> 00:01:16,978
We have seen that the simplified Newton method is a reasonable approach to the iterative

16
00:01:17,145 --> 00:01:20,025
solution of these increments equations.

17
00:01:20,192 --> 00:01:26,472
However if it takes several iterations to solve them this can still be very expensive.

18
00:01:26,639 --> 00:01:31,719
And now in the current section we are looking at so called semi implicit Runge Kutta methods

19
00:01:31,886 --> 00:01:39,406
which try to avoid many iterations required for the solution of the increment equations.

20
00:01:39,573 --> 00:01:41,053
So what is the idea?

21
00:01:41,220 --> 00:01:42,820
The idea is fairly natural.

22
00:01:42,987 --> 00:01:49,227
The idea is to use only a fixed small number of iterations and there is a natural motivation

23
00:01:49,394 --> 00:01:50,554
underlying this idea.

24
00:01:50,721 --> 00:01:56,041
Remember that it is absolutely pointless to solve the increment equations up to machine

25
00:01:56,208 --> 00:02:02,048
precision because the single step method will only provide an approximate solution anyway.

26
00:02:02,215 --> 00:02:06,595
A solution states which have a discretization error.

27
00:02:06,762 --> 00:02:14,662
So the iterative solver need only solve down to the level of that discretization error.

28
00:02:14,829 --> 00:02:19,909
The second reasoning is that we already have a good initial guess for the iteration.

29
00:02:20,076 --> 00:02:25,316
When we look at the increment equations in the stage form then we can use 0 at initial

30
00:02:25,483 --> 00:02:30,763
guess and that initial guess will be the closer to the final solution the smaller the time

31
00:02:30,930 --> 00:02:32,290
step.

32
00:02:32,457 --> 00:02:34,057
There is a nice balancing.

33
00:02:34,224 --> 00:02:38,064
Small time step means small discretization error on the one hand.

34
00:02:38,231 --> 00:02:45,071
On the other hand it also means an excellent initial guess for the solution of the nonlinear

35
00:02:45,238 --> 00:02:46,518
system.

36
00:02:46,685 --> 00:02:54,165
This justifies the hope that maybe a small fixed number or in the extreme case a single

37
00:02:54,332 --> 00:02:56,012
Newton step will suffice.

38
00:02:56,179 --> 00:03:01,819
So let's first try this idea in a concrete setting and carry out a numerical experiment.

39
00:03:01,986 --> 00:03:05,406
We apply the implicit Euler method to an autonomous ODE.

40
00:03:05,573 --> 00:03:12,413
So as all of you know this yields this recursion for advancing from Tk to Tk plus 1.

41
00:03:12,580 --> 00:03:19,260
And now let's figure out how the method will be modified if we rely on a single Newton

42
00:03:19,427 --> 00:03:21,987
step to solve this nonlinear equation.

43
00:03:22,154 --> 00:03:27,894
To that end we convert it into an equivalent zero finding problem for a function capital

44
00:03:28,061 --> 00:03:34,701
F. So F of Yk plus 1 here is so we bring everything to one side.

45
00:03:34,868 --> 00:03:45,268
Let's do it like this minus Yk minus Hf of Yk plus 1 equal to 0 and solve this for Yk

46
00:03:45,435 --> 00:03:46,435
plus 1.

47
00:03:46,602 --> 00:03:52,842
Now let's apply a single Newton step with initial guess Yk which is the natural initial

48
00:03:53,009 --> 00:03:54,689
guess for the next state.

49
00:03:54,856 --> 00:03:56,316
So what do we get?

50
00:03:56,483 --> 00:03:59,183
So you know the formula for the Newton iteration.

51
00:03:59,350 --> 00:04:01,910
So this is just one Newton step.

52
00:04:02,077 --> 00:04:09,917
So starting from Yk and then it will be the linear system involving the Jacobian of capital

53
00:04:10,084 --> 00:04:17,504
F. Apply to F of Yk and now it's easy to compute the Jacobian of F with respect to its

54
00:04:17,671 --> 00:04:18,531
argument.

55
00:04:18,698 --> 00:04:24,018
This is an n by n matrix if you're in n dimensional phase space.

56
00:04:24,185 --> 00:04:28,605
So let's differentiate F with respect to Yk plus 1.

57
00:04:28,772 --> 00:04:38,752
This is minus HdFYk plus 1 inverse and now F of Yk becomes simple.

58
00:04:38,919 --> 00:04:46,439
It's minus Hf of Yk plus 1 so take into account everything by putting a plus here.

59
00:04:46,606 --> 00:04:49,286
This is the new recursion and this is a known method.

60
00:04:49,453 --> 00:04:52,573
This is the so called semi implicit Euler method.

61
00:04:52,740 --> 00:04:55,620
Well it was easy to derive it but will it perform?

62
00:04:55,787 --> 00:04:59,327
Let's study this in a numerical experiment.

63
00:04:59,494 --> 00:05:04,674
So we apply the semi implicit Euler method to an initial value problem for the logistic

64
00:05:04,841 --> 00:05:11,361
differential equation and this is a moderately stiff initial value problem because the lambda

65
00:05:11,528 --> 00:05:12,928
is said to be 5.

66
00:05:13,095 --> 00:05:17,735
We run it with uniform time step for many different time steps that we obtain just by

67
00:05:17,902 --> 00:05:21,462
halving the time step size repeatedly.

68
00:05:21,629 --> 00:05:27,689
And we measure the error in the points of the temporal mesh and take the maximum.

69
00:05:27,856 --> 00:05:33,876
And this is what is plotted here for the implicit Euler in blue and for semi implicit Euler

70
00:05:34,043 --> 00:05:35,603
in red.

71
00:05:35,770 --> 00:05:37,410
Double logarithmic plot.

72
00:05:37,577 --> 00:05:46,177
We clearly see algebraic convergence but the great news is that both methods enjoy algebraic

73
00:05:46,344 --> 00:05:48,584
convergence with rate 1.

74
00:05:48,751 --> 00:05:57,011
So fixing the number of Newton steps to only 1 has not led to any deterioration.

75
00:05:57,178 --> 00:06:02,718
Neither in terms of rate nor as you can see in terms of actual accuracy.

76
00:06:02,885 --> 00:06:11,365
So we still achieve the same order and accuracy with the semi implicit Euler method as with

77
00:06:11,532 --> 00:06:12,792
the fully implicit.

78
00:06:12,959 --> 00:06:14,539
Let me remark this here.

79
00:06:14,706 --> 00:06:20,666
Well it seems that our initial reasoning was not really off the mark.

80
00:06:20,833 --> 00:06:21,833
Was a valid argument.

81
00:06:22,000 --> 00:06:25,120
Of course implicit Euler is a really low order method.

82
00:06:25,287 --> 00:06:26,287
Order 1.

83
00:06:26,454 --> 00:06:29,754
So let's see whether the idea carries over to methods of high order.

84
00:06:29,921 --> 00:06:35,801
So we look at the semi implicit midpoint method which we obtain from the fully implicit midpoint

85
00:06:35,968 --> 00:06:40,468
method by solving the increment equation with a single Newton step.

86
00:06:40,635 --> 00:06:44,915
Again same model problem and this is the result.

87
00:06:45,082 --> 00:06:49,242
Uniform time steps and the same error.

88
00:06:49,409 --> 00:06:53,949
Again algebraic convergence is really conspicuous in this doubly logarithmic plot.

89
00:06:54,116 --> 00:07:01,576
We see that again the algebraic rate of convergence is preserved.

90
00:07:01,743 --> 00:07:07,303
So using a single Newton step apparently has not destroyed the second order of the method.

91
00:07:07,470 --> 00:07:10,590
We incur a slight loss of accuracy but not severely.

92
00:07:10,757 --> 00:07:17,197
So let's summarize that the semi implicit midpoint method still enjoys order 2 of the

93
00:07:17,364 --> 00:07:19,604
fully implicit midpoint method.

94
00:07:19,771 --> 00:07:21,451
Well this is encouraging.

95
00:07:21,618 --> 00:07:27,338
Could it be that this simple approach remains valid for all implicit Runge Kutta methods

96
00:07:27,505 --> 00:07:28,685
of arbitrary order.

97
00:07:28,852 --> 00:07:35,112
So can it be that this the use of a single Newton step which is equivalent to the plain

98
00:07:35,279 --> 00:07:41,459
and simple linearization of the increment equation also works for high order methods

99
00:07:41,626 --> 00:07:43,566
without destroying their order.

100
00:07:43,733 --> 00:07:44,973
Let's try.

101
00:07:45,140 --> 00:07:49,620
So let's try to use linearized increment equations.

102
00:07:49,787 --> 00:07:54,687
And this is more or less equivalent to the use of a single Newton step.

103
00:07:54,854 --> 00:08:01,554
So what we do is we linearize the increment equations by expanding f in well a Taylor

104
00:08:01,721 --> 00:08:05,121
series and truncating after the linear term.

105
00:08:05,288 --> 00:08:08,008
We expand around y naught.

106
00:08:08,175 --> 00:08:14,715
So now I consider the first step from y naught to y1 and then this is the Taylor expansion.

107
00:08:14,882 --> 00:08:20,222
This is the sort of perturbation and this again shows up here.

108
00:08:20,389 --> 00:08:26,949
And all say terms that are quadratic in the step size h are dropped.

109
00:08:27,116 --> 00:08:29,696
Now let's see where this idea can work.

110
00:08:29,863 --> 00:08:31,003
Some observations.

111
00:08:31,170 --> 00:08:37,090
In the case of a linear ODE or even affine linear ODE the linearization of course has

112
00:08:37,257 --> 00:08:44,117
no effect because in this case the increment equations are linear anyway.

113
00:08:44,284 --> 00:08:51,264
Since the stability analysis was based on studying the single step method for linear

114
00:08:51,431 --> 00:08:57,151
ODE's all what we have learned through this linear model problem analysis remains valid.

115
00:08:57,318 --> 00:09:00,758
In particular the stability functions will not change.

116
00:09:00,925 --> 00:09:07,525
So if we have an unconditionally stable implicit Runge Kutta single step method then what we

117
00:09:07,692 --> 00:09:12,092
obtain after linearization will still enjoy that property.

118
00:09:12,259 --> 00:09:18,459
Another crucial question as I have already pointed out is the question whether the order

119
00:09:18,626 --> 00:09:20,706
of the method is affected.

120
00:09:20,873 --> 00:09:23,253
Let's look at even higher order case.

121
00:09:23,420 --> 00:09:31,260
The semi implicit two stage Raudau method again applied to our standard model problem

122
00:09:31,427 --> 00:09:32,407
for this section.

123
00:09:32,574 --> 00:09:39,174
The two stage Raudau method is an implicit Runge Kutta single step method of order three.

124
00:09:39,341 --> 00:09:46,541
So we do the linearization, studied the error and its convergence as we send the step size,

125
00:09:46,708 --> 00:09:50,828
the uniform step size to zero and now we are disappointed.

126
00:09:50,995 --> 00:09:57,515
What we clearly see is that the method obtained by linearization still shows a nice algebraic

127
00:09:57,682 --> 00:10:03,442
convergence but of lower order than the original Runge Kutta single step method.

128
00:10:03,609 --> 00:10:09,409
Two stage Raudau order three, naive semi implicit version only order two.

129
00:10:09,576 --> 00:10:14,176
So we have to concede a loss of order due to linearization.

130
00:10:14,343 --> 00:10:16,463
Well should we give up the idea?

131
00:10:16,630 --> 00:10:17,130
No.

132
00:10:17,297 --> 00:10:23,277
We rather use this idea of relying on a single Newton step as an inspiration.

133
00:10:23,444 --> 00:10:24,444
Look to the left.

134
00:10:24,611 --> 00:10:31,931
The linearized increment equations still contain the usual Runge Kutta coefficients and now

135
00:10:32,098 --> 00:10:37,978
researchers in numerical analysis said to themselves let's pursue an approach similar

136
00:10:38,145 --> 00:10:43,985
to the say the bootstrapping construction of Runge Kutta methods or the collocation

137
00:10:44,152 --> 00:10:47,052
construction of implicit Runge Kutta methods.

138
00:10:47,219 --> 00:10:53,999
Those techniques gave us formulas that had a special structure with particular parameters

139
00:10:54,166 --> 00:11:00,326
but then we just started from those formulas and used them as a template to build more

140
00:11:00,493 --> 00:11:08,373
methods and now we do the same here and this leads us to the so called class of Rosenbrock

141
00:11:08,540 --> 00:11:09,580
Banach methods.

142
00:11:09,747 --> 00:11:15,187
So for a general autonomous ODE we write down these formulas.

143
00:11:15,354 --> 00:11:17,554
They are inspired by linearization.

144
00:11:17,721 --> 00:11:19,881
You recognize the Jacobian here.

145
00:11:20,048 --> 00:11:25,048
These formulas are closely related to the linearized increment equations but now we

146
00:11:25,215 --> 00:11:31,615
have introduced just more parameters and now what has been done is people started from

147
00:11:31,782 --> 00:11:39,702
these formulas and then they determine the parameters or coefficients the aij's or the

148
00:11:39,869 --> 00:11:45,789
dij's to achieve a prescribed order and to ensure some stability properties.

149
00:11:45,956 --> 00:11:49,676
So this is a stability, l stability.

150
00:11:49,843 --> 00:11:55,823
This is not easy because to achieve a prescribed order you have to solve so called order conditions

151
00:11:55,990 --> 00:12:01,010
which are large nonlinear systems for these coefficients.

152
00:12:01,177 --> 00:12:09,017
However it's feasible and this gives rise to a large family of methods that still feature

153
00:12:09,184 --> 00:12:13,664
nice stability properties so they are suitable for stiff problems but we no longer have to

154
00:12:13,831 --> 00:12:16,311
solve nonlinear systems in each step.

155
00:12:16,478 --> 00:12:21,518
All we have to do is to solve these equations and what does this boil down to?

156
00:12:21,685 --> 00:12:29,125
This is an, well an n by n matrix so we have to solve an n by n linear system of equations

157
00:12:29,292 --> 00:12:37,252
in each step to determine the increment vectors and you can imagine that the relative straightforward

158
00:12:37,419 --> 00:12:43,499
implementation of a single step combined with good stability properties makes this class

159
00:12:43,666 --> 00:12:46,666
very attractive for practical applications.

160
00:12:46,833 --> 00:12:53,113
So Rosenbrock Wanner method, abbreviated as Rowe methods are nowadays widely used in practical

161
00:12:53,280 --> 00:12:58,000
application for the simulation of time dependent phenomena.