Combine linear maps in interval-matrix discretization

[schillic]

Instead of taking the sum of three LinearMap instances all with the set U0, we can add up the matrices and take the LinearMap in the end. This is both more efficient and more precise (because we need to overapproximate only once).

Reachability.jl/src/ReachSets/discretize.jl

Lines 803 to 825 in e792eb0

	A² = A * A
	Iδ = IntervalMatrix(Diagonal(fill(IntervalMatrices.Interval(δ), size(A, 1))))
	IδW = Iδ + 1/2 * δ^2 * A + 1/6 * δ^3 * A²
	linear_maps[1] = LinearMap(IδW, U0)
	E = _expm_remainder(A, δ, order; n=n)
	linear_maps[2] = LinearMap(E*δ, U0)
	zero_interval = IntervalMatrices.Interval(zero(N), zero(N))
	if order > 2
	# i = 2
	αᵢ₊₁ = 6 # factorial of (i+1)
	Aⁱ = A²
	δⁱ⁺¹ = δ^3
	M_sum = IntervalMatrix(fill(zero_interval, size(A)))
	@inbounds for i in 3:order
	αᵢ₊₁ *= i+1
	δⁱ⁺¹ *= δ
	Aⁱ *= A
	M_sum += (δⁱ⁺¹/αᵢ₊₁) * Aⁱ
	end
	linear_maps[3] = LinearMap(M_sum, U0)
	end