BasicVariable.m 7.46 KB
Newer Older
1
classdef BasicVariable < handle
2
3
4
5
6
	%BasicVariable Class to handle quantities with a priori known derivatives
	% This class is used to simplify the usage of quantities for which a lot of computations must be
	% done on its derivatives. So a BasicVariable consists of a basic function and its derivatives,
	% which are precomputed and safed. So computations can be performed without recalculation of the
	% derivatives.
7
8
	
	properties ( SetAccess = protected )
9
10
		derivatives cell;
		fun (1,1) quantity.Discrete;
11
12
	end
	
13
14
15
16
17
18
	properties ( Dependent )
		highestDerivative;  % order of the highest computed derivative
		T;					% the transition time T
		dt;					% step size
	end
	
19
	methods
20
21
22
23
		function obj = BasicVariable(fun, derivatives)
			obj.fun = fun;
			obj.derivatives = derivatives;
		end
24
25
26
27
		function h = get.highestDerivative(obj)
			h = numel(obj(1).derivatives);
		end
		function T = get.T(obj)
28
			T = obj.domain.upper();
29
30
		end
		function dt = get.dt(obj)
31
32
			if isa(obj.domain, "quantity.EquidistantDomain")
				dt = obj.domain.stepSize;
33
			else
34
				delta_t = diff(obj.domain.grid);
35
36
37
38
				assert( all( diff( delta_t ) < 1e-15 ), "Grid is not equidistant spaced" );
				dt = delta_t(1);
			end
		end
39
40
41
	end
	
	methods ( Access = public)
42
43
44
45
46
47
48
49
		function D = diff(obj, domain, n)
			% DIFF get the k-th derivative 
			%	D = diff(obj, domain, n) tries to find the n-th derivative of obj.fun in the saved
			%	derivatives. If it does not exist, it tries to compute it with the evaluateFunction.
			%	If n is a non integer number, the Riemann-Liouville fractional derivative is
			%	computed. For this, a formula is used which requires the next integer + 1 derivative
			%	of the function to avoid the singularity under the integral. Actually, the domain is
			%	clear, but it is used hear to be consistent with the quantity.Discrete/diff.
50
51
			arguments
				obj
52
53
				domain (1,1) quantity.Domain;
				n (1,1) double = 0;
54
			end
55
			
56
			for l = 1:numel(obj)
57
58
59
60
61
				
				assert( domain.isequal( obj(l).domain  ),...
					"The derivative wrt. this domain is not possible.")
				
				if n == 0
62
					D(l) = obj(l).fun;
63
64
65
66
67
68
					
				elseif ~numeric.near(round(n), n)
					% fractional derivative:
					D(l) = obj(l).fDiff(n);
					
				elseif n > length(obj(l).derivatives)
69
70
71
					% For some functions their maybe a possibility to compute new derivatives on the
					% fly. For this, the evaluateFunction(k) can be used, where k is the order of
					% the derivative.
72
					D(l) = obj(l).evaluateFunction(obj.domain, n);
73
					obj(l).derivatives{end+1} = D(l);
74
				else
75
76
					% use the pre-computed derivatives
					D(l) = obj(l).derivatives{n};
77
78
79
80
81
82
83
84
85
86
87
88
89
				end
			end
			D = reshape(D, size(obj));
		end
		
		function D = diffs(obj, k)
			% DIFFS compute multiple derivatives
			%	D = diffs(obj, k) computes the k-th derivatives of OBJ. At this, k can be a vector
			%	of derivative orders.
			arguments
				obj,
				k (:,1) double {mustBeInteger, mustBeNonnegative} = 0:numel(obj(1).derivatives);
			end
90
			D_ = arrayfun( @(i) obj.diff(obj.domain, i), k, 'UniformOutput', false);
91
			
92
93
			for i = 1:numel(k)
				D(i,:) = D_{i};
94
95
96
			end
		end
		
97
98

		
99
100
101
102
103
104
105
106
107
108
109
110
111
		function h = productRule(obj, phi, n)
			% LEIBNIZ compute the k-derivative of the product of the basic variable and phi
			%	h = productRule(obj, phi, n) computes the k-th derivative of obj.fun*phi, i.e.,
			%		h = d_t^n ( obj.fun(t) * phi ) 
			%	using the product rule. This is usefull if the derivatives of the basic variable obj
			%	are known exactly and the derivatives of phi, but not the multiplication
			%	obj.diff(0)*phi.
			t = phi.domain;
			h = quantity.Discrete.zeros( size(obj), t );
			
			if numeric.near(round(n), n)
				% for integer order derivatives
				for k = 0:n
112
					h = h + misc.binomial(n, k) * obj.diff(t, n-k) * phi.diff(t, k);
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
				end
			else
				% for fractional order derivatives: [see Podlubny]
				%         p                     \~~ oo  / p \     (k)       p-k
				%        D  (phi(t) f(t)) =  >      |   |    | phi   (t)   D   f(t)
				%      a  t                     /__ k=0 \ k /               a  t				
				h = misc.binomial( n, 0 ) * phi * obj.fDiff(n);
				h0 = h;
				k = 1;
				% do the iteration on the loop until the difference is very small.
				while max( abs( h0 ) ) / max( abs( h ) ) > 1e-9 
					h0 = misc.binomial(n,k) * phi.diff(t, k) * obj.fDiff( n - k);
					h = h + h0;
					k = k+1;
				end					
			end
129
130
131
		end
		
		
132
		function d = domain(obj)
133
			d = obj(1).fun(1).domain;
134
		end
135
136
		
		function plot(obj, varargin)
137
			F = quantity.Discrete(obj);			
138
139
			F.plot(varargin{:});
		end
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
		
		function F = quantity.Discrete(obj)
			F (1,:) = [obj.fun];
			
			for i = 1:max([obj.highestDerivative])
							
				for j = 1:numel(obj)
					F_i(j) = obj(j).derivatives{i};
				end
				
				F(i+1,:) = reshape(F_i, size(obj));
			end
		end		
		
		function c = mtimes(a, b)
			arguments
				a double;
				b signals.BasicVariable
			end
			
			assert( size(a,2) == size(b,1), "dimensions of the terms for multiplication do not match")
			F = quantity.Discrete(b);
			
			% do the multiplication for each derivative:
			for i = 1:max([b.highestDerivative])+1
				tmp(i,:) = a * shiftdim(F(i,:));
			end
			
			% restore the result as signals.BasicVariable
			for k = 1:numel(b)
				for l = 2:size(tmp,1)
					derivs_{l-1} = tmp(l,k);
				end
				c(k) = signals.BasicVariable( tmp(1,k), derivs_);
			end
			c = reshape(c, size(b));
		end
			
	end % methods (Access = public)
	
	methods (Access = protected)
181
		function v = evaluateFunction(obj, domain, k)
182
183
			error('conI:signals:BasicVariable:derivative', ['Requested derivative ' num2str(k) ' is not available']);
		end
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
		
		function D = fDiff(obj, p)
			% FDIFF compute the p-fractional derivative
			%	D = fDiff(obj, p) computes the fractional derivative of order p. To be
			%	concrete the Riemann-Liouville fractional derivative with initial point t0 = lower
			%	end of the considered time domain is computed.
			%	Because the basisc variable must be C^infinity function, we use the form
			%	of the fractional derivative described in eq. (2.80) in the book [Podlubny: Fractional
			%	differential equations]. Then, integration by parts is applied to eliminate the
			%	singularity:
			%	
			
			n = ceil(p);
			tDomain = obj(1).domain;
			t0 = tDomain.lower;
			t = Discrete( tDomain );
			tauDomain = tDomain.rename( tDomain.name + "_Hat" );
			tau = Discrete( tauDomain );
			if p > 0
				
				% TODO Try to implement the polynomial part as function or as symbolic
				D = quantity.Discrete.zeros(size(obj), tDomain);
				for k = 0:n
					f_dk = obj.diff(tDomain, k);
					
					if f_dk.at(t0) ~= 0
						D = D + f_dk.at(t0) * (t - t0).^(k-p) / gamma( k-p + 1);
					end
					
				end
				assert( ~ any( isinf( obj.diff(tDomain, n+1).on() ) ), "conI:BasicVariable:fDiff", ...
					"The function has a singularity! Hence, the fractional derivative can not be computed")
				D = D + cumInt( ( t - tau ).^(n-p) * subs(obj.diff(tDomain, n+1), tDomain.name, tauDomain.name), ...
				tauDomain, tauDomain.lower, tDomain.name) / gamma( n+1-p );
				
			elseif p < 0 && p > -1
				alpha = -p;
				D = (t - t0).^(alpha) / gamma( alpha + 1) * obj.diff(tDomain, 0).at(t0) + ...
					cumInt( (t-tau).^alpha * subs(obj.diff(tDomain, 1), tDomain.name, tauDomain.name), ...
					tauDomain, tauDomain.lower, tDomain.name) / gamma( alpha + 1 );

			elseif p <= -1
				alpha = -p;
				D = cumInt( (t-tau).^(alpha-1) * subs( obj.diff(tDomain, 0), tDomain.name, tauDomain.name), ...
					tauDomain, tauDomain.lower, tDomain.name) / gamma( alpha );
				
			elseif p == 0
				D = [obj.fun()];	
			else
				error("This should not happen")
			end
		end		
		
237
238
	end	
	
239
end