Skip to content
Snippets Groups Projects
Commit 2589646e authored by Gfrerer, Michael Helmut's avatar Gfrerer, Michael Helmut
Browse files

hyper dual numbers added

parent 97a9aa6c
Branches
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Showing
with 567 additions and 0 deletions
matlab/hyper_dual_numbers/@hyper_dual/horzcat.m 0 → 100644
+ 4
− 0
View file @ 2589646e
function ret = horzcat(varargin)
obj = cat( 1, varargin{:} );
ret = hyper_dual([obj.a],[obj.b],[obj.c],[obj.d]);
end
matlab/hyper_dual_numbers/@hyper_dual/hyper_dual.m 0 → 100644
+ 439
− 0
View file @ 2589646e
classdef hyper_dual
%HYPER_DUAL
properties
a % the number
b
c
d
end
methods
function obj = hyper_dual(A,B,C,D)
if nargin == 1
obj.a = zeros(size(A));
obj.b = zeros(size(A));
obj.c = zeros(size(A));
obj.d = zeros(size(A));
obj.a = A;
elseif nargin == 4
if isequal(size(A),size(B)) && isequal(size(B),size(C)) && isequal(size(C),size(D))
obj.a = A;
obj.b = B;
obj.c = C;
obj.d = D;
else
error('wrong input')
end
else
error('wrong number of input parameter in hyper_dual')
end
end
%% double
function ret = double(obj)
ret = obj;
ret.a = double(obj.a);
ret.b = double(obj.b);
ret.c = double(obj.c);
ret.d = double(obj.d);
end
%% real
function ret = real(obj)
ret = obj.a;
end
%% imag
function ret = imag(obj)
ret = abs(imag(obj.b))+abs(imag(obj.c))+abs(imag(obj.d));
end
%% reshape
function ret = reshape(obj,varargin)
ret = hyper_dual(...
reshape(obj.a,varargin{:}),...
reshape(obj.b,varargin{:}),...
reshape(obj.c,varargin{:}),...
reshape(obj.d,varargin{:}));
end
function ret = transpose(obj)
ret = hyper_dual(obj.a.',obj.b.',obj.c.',obj.d.');
end
%% addition
function ret = plus(obj,A)
if isa(obj , 'hyper_dual')
ret = obj;
if isa(A , 'hyper_dual')
ret.a = bsxfun(@plus,ret.a , A.a); % if A is hyper-dual
ret.b = bsxfun(@plus,ret.b , A.b);
ret.c = bsxfun(@plus,ret.c , A.c);
ret.d = bsxfun(@plus,ret.d , A.d);
else
ret.a = bsxfun(@plus,ret.a , A);
end
else
ret = A + obj;
end
end
%% subtraction
function ret = minus(obj,A)
if isa(obj , 'hyper_dual')
ret = obj;
if isa(A , 'hyper_dual')
ret.a = bsxfun(@minus,ret.a , A.a); % if A is hyper-dual
ret.b = bsxfun(@minus,ret.b , A.b);
ret.c = bsxfun(@minus,ret.c , A.c);
ret.d = bsxfun(@minus,ret.d , A.d);
else
ret.a = bsxfun(@minus,ret.a , A);
end
else
ret = - A + obj ;
end
end
%% minus
function ret = uminus(obj)
ret = obj;
f = -obj.a;
df = -1+0*obj.a;
ddf= 0*obj.a;
ret.a = f; % if A is hyper-dual
ret.b = obj.b .* df;
ret.c = obj.c .* df;
ret.d = obj.d .* df + obj.b .* obj.c .* ddf;
end
%% Multiplication
function ret = times(obj,A)
if isa(obj , 'hyper_dual')
ret = obj;
if isa(A , 'hyper_dual')
ret.a = bsxfun(@times,obj.a , A.a); % if A is hyper-dual
ret.b = bsxfun(@times,obj.a , A.b) + bsxfun(@times,obj.b , A.a);
ret.c = bsxfun(@times,obj.a , A.c) + bsxfun(@times,obj.c , A.a);
ret.d = bsxfun(@times,obj.a , A.d) + bsxfun(@times,obj.b , A.c) + bsxfun(@times,obj.c , A.b) + bsxfun(@times,obj.d , A.a);
else
ret = obj;
ret.a = bsxfun(@times,obj.a , A); % if A is a real number
ret.b = bsxfun(@times,obj.b , A);
ret.c = bsxfun(@times,obj.c , A);
ret.d = bsxfun(@times,obj.d , A);
end
else
ret = A .* obj;
end
end
%% Matrix Multiplication
function ret = reini(obj)
for i = 1:size(obj,1)
for j = 1:size(obj,2)
s1(i,j) = size(obj(i,j).a,1);
s2(i,j) = size(obj(i,j).a,2);
end
end
i1 = 0;
for i = 1:size(obj,1)
i2 = 0;
for j = 1:size(obj,2)
for k = 1:s1(i,j)
for l = 1:s2(i,j)
a(i1+k,i2+l) = obj(i,j).a(k,l);
b(i1+k,i2+l) = obj(i,j).b(k,l);
c(i1+k,i2+l) = obj(i,j).c(k,l);
d(i1+k,i2+l) = obj(i,j).d(k,l);
end
end
i2 = i2 + s2(i,j);
end
i1 = i1 + s1(i,j);
end
ret = hyper_dual(a,b,c,d);
end
%%
function ret = mtimes(obj,A)
if isa(obj , 'hyper_dual')
if numel(obj) > 1
obj = reini(obj);
end
else
obj = hyper_dual(obj);
end
sO1 = size(obj.a,1);
sO2 = size(obj.a,2);
if isa(A , 'hyper_dual')
if numel(A) > 1
A = reini(A);
end
else
A = hyper_dual(A);
end
sA1 = size(A.a,1);
sA2 = size(A.a,2);
if sO2 == sA1
aa = zeros(sO1,sA2);
bb = zeros(sO1,sA2);
cc = zeros(sO1,sA2);
dd = zeros(sO1,sA2);
for i = 1:sO1
for j = 1:sA2
for k = 1:sO2
r = my_subs(obj,i,k) .* my_subs(A,k,j);
aa(i,j) = aa(i,j) + r.a;
bb(i,j) = bb(i,j) + r.b;
cc(i,j) = cc(i,j) + r.c;
dd(i,j) = dd(i,j) + r.d;
end
end
end
ret = hyper_dual(aa,bb,cc,dd);
else
error('dimension mismatch')
end
end
%% DIVIDE
% rdivide ./
% obj is nominator
% A is denominator
function ret = rdivide(obj,A)
if isa(obj , 'hyper_dual')
ret = obj;
if isa(A , 'hyper_dual')
df = bsxfun(@rdivide,-1, A.a.^2);
ddf= bsxfun(@rdivide,2, A.a.^3);
ret.a = bsxfun(@rdivide,1, A.a); % if A is hyper-dual
ret.b = A.b .* df;
ret.c = A.c .* df;
ret.d = A.d .* df + A.b .* A.c .* ddf;
ret = obj .* ret;
else
ret = obj .* (1 ./ A);
end
else
ret = hyper_dual(obj) ./ A;
end
end
%% Equality
function ret = eq(obj,a)
if isa(obj , 'hyper_dual')
if isa(a , 'hyper_dual')
% if a is hyper-dual
ret = (obj.a == a.a) ;
else
ret = ( obj.a == a ); % if a is a real number
end
else
ret = ( obj == a.a );
end
end
%% sin
function ret = sin(obj)
ret = obj;
ret.a = sin( obj.a );
ret.b = obj.b .* cos(obj.a);
ret.c = obj.c .* cos(obj.a);
ret.d = obj.d .* cos(obj.a) - obj.b .* obj.c .* sin( obj.a );
end
%% cos
function ret = cos(obj)
ret = obj;
ret.a = cos( obj.a );
ret.b = -obj.b .* sin(obj.a);
ret.c = -obj.c .* sin(obj.a);
ret.d = -obj.d .* sin(obj.a) - obj.b .* obj.c .* cos( obj.a );
end
%% exp
function ret = exp(obj)
ret = obj;
ret.a = exp( obj.a );
ret.b = obj.b .* exp(obj.a);
ret.c = obj.c .* exp(obj.a);
ret.d = obj.d .* exp(obj.a) + obj.b .* obj.c .* exp(obj.a);
end
% power
function ret = power(obj,B)
if isa(obj , 'hyper_dual')
if isa(B , 'hyper_dual')
% if B is hyper-dual
error('not implemented')
else
% if B is a real number
ret = obj;
help = B .* obj.a .^ (B-1);
ret.a = obj.a .^ B;
ret.b = obj.b .* help;
ret.c = obj.c .* help;
ret.d = obj.d .* help + obj.b .* obj.c .* B .* (B-1) .* obj.a .^ (B-2);
end
else
ret = B; % B is hyper-dual
ret.a = obj .^ B.a;
ret.b = B.b .* obj .^ B.a .* log(obj) ;
ret.c = B.c .* obj .^ B.a .* log(obj) ;
ret.d = B.d .* obj .^ B.a .* log(obj) + B.b .* B.c .* obj .^ B.a .* log(obj)^2;
end
end
%% sqrt
function ret = sqrt(obj)
ret = obj;
Df = 1 ./ ( 2 * sqrt(obj.a) );
ret.a = sqrt(obj.a);
ret.b = obj.b .* Df;
ret.c = obj.c .* Df;
ret.d = obj.d .* Df - obj.b .* obj.c .* Df ./ ( 2 * obj.a );
end
%% sinh
function ret = sinh(obj)
ret = obj;
Df = cosh(obj.a) ;
ret.a = sinh(obj.a);
ret.b = obj.b .* Df;
ret.c = obj.c .* Df;
ret.d = obj.d .* Df + obj.b .* obj.c .* sinh(obj.a);
end
%% cosh
function ret = cosh(obj)
ret = obj;
Df = sinh(obj.a) ;
ret.a = cosh(obj.a);
ret.b = obj.b .* Df;
ret.c = obj.c .* Df;
ret.d = obj.d .* Df + obj.b .* obj.c .* cosh(obj.a);
end
%% abs
function ret = abs(obj)
ret = obj;
ret.a = abs(obj.a);
end
%% <
function ret = lt(obj,B)
if isa(obj , 'hyper_dual')
if isa(B , 'hyper_dual')
% obj and B is hyper-dual
ret = obj.a < B.a;
else
% if B is a real number
ret = obj.a < B;
end
else
% B is hyper-dual
ret = obj < B.a;
end
end
%% <=
function ret = le(obj,B)
if isa(obj , 'hyper_dual')
if isa(B , 'hyper_dual')
% obj and B is hyper-dual
ret = obj.a <= B.a;
else
% if B is a real number
ret = obj.a <= B;
end
else
% B is hyper-dual
ret = obj <= B.a;
end
end
%% >
function ret = gt(obj,B)
if isa(obj , 'hyper_dual')
if isa(B , 'hyper_dual')
% obj and B is hyper-dual
ret = obj.a > B.a;
else
% if B is a real number
ret = obj.a > B;
end
else
% B is hyper-dual
ret = obj > B.a;
end
end
%% >=
function ret = ge(obj,B)
if isa(obj , 'hyper_dual')
if isa(B , 'hyper_dual')
% obj and B is hyper-dual
ret = obj.a >= B.a;
else
% if B is a real number
ret = obj.a >= B;
end
else
% B is hyper-dual
ret = obj >= B.a;
end
end
%% mod
function ret = mod(obj,B)
if isa(obj , 'hyper_dual')
if isa(B , 'hyper_dual')
% obj and B is hyper-dual
error('not implemented')
ret = mod(obj.a , B.a);
else
% if B is a real number
ret = obj;
ret.a = mod(obj.a , B);
end
else
% B is hyper-dual
error('not implemented')
ret = mod(obj , B.a);
end
end
%% sign
function ret = sign(obj)
ret.a = sign(obj.a);
end
% function ret = subsasgn(obj,s,a)
% if length(s) == 1
% if s.type == '()'
% error('not implemented')
% else
% ret = builtin('subsasgn', obj, s,a);
% end
% elseif length(s) == 2
% if strcmp(s(1).type,'.') && strcmp(s(2).type,'()')
% ret = builtin('subsasgn', obj, s,a);
% else
% error('not implemented')
% end
% else
% error('not implemented')
% end
% end
%function C = horzcat(obj,B)
%end
%%
%%
function n = numArgumentsFromSubscript(obj, s, ic)
n = builtin('numArgumentsFromSubscript', obj, s, ic);
end
end
methods (Static)
%% zeros
function ret = zeros(numX,numY)
ret = hyper_dual( zeros(numX,numY) );
end
end
methods (Hidden)
%% zeros
function ret = zerosLike(obj,numX,numY)
if nargin == 2
ret = hyper_dual( zeros(numX) );
else
ret = hyper_dual( zeros(numX,numY) );
end
end
%% my_subs()
function ret = my_subs(obj,i,j)
ret = hyper_dual(obj.a(i,j),obj.b(i,j),obj.c(i,j),obj.d(i,j));
end
end
end
matlab/hyper_dual_numbers/@hyper_dual/subsref.m 0 → 100644
+ 27
− 0
View file @ 2589646e
function ret = subsref(obj,s)
ret = obj;
if length(s) == 1
if strcmp(s.type,'()')
if length(s.subs) == 1
ret.a = obj.a(s.subs{1});
ret.b = obj.b(s.subs{1});
ret.c = obj.c(s.subs{1});
ret.d = obj.d(s.subs{1});
elseif length(s.subs) == 2
ret.a = obj.a(s.subs{1},s.subs{2});
ret.b = obj.b(s.subs{1},s.subs{2});
ret.c = obj.c(s.subs{1},s.subs{2});
ret.d = obj.d(s.subs{1},s.subs{2});
else
ret.a = builtin('subsref',obj.a,s);
ret.b = builtin('subsref',obj.b,s);
ret.c = builtin('subsref',obj.c,s);
ret.d = builtin('subsref',obj.d,s);
end
else
ret = builtin('subsref', obj, s);
end
else
ret = builtin('subsref', obj, s);
end
end
\ No newline at end of file
matlab/hyper_dual_numbers/@hyper_dual/vertcat.m 0 → 100644
+ 16
− 0
View file @ 2589646e
function ret = vertcat(varargin)
obj = cat( 1, varargin{:} );
ret = hyper_dual(cat(1,obj.a),cat(1,obj.b),cat(1,obj.c),cat(1,obj.d));
% if isa(A , 'hyper_dual')
% if isa(B , 'hyper_dual')
% % obj and B is hyper-dual
% ret = hyper_dual([A.a;B.a],[A.b;B.b],[A.c;B.c],[A.d;B.d]);
% else
% % if B is a real number
% ret = hyper_dual([A.a;B],[A.b;0*B],[A.c;0*B],[A.d;0*B]);
% end
% else
% % B is hyper-dual
% ret = hyper_dual([A;B.a],[0*A;B.b],[0*A;B.c],[0*A;B.d]);
% end
end
matlab/hyper_dual_numbers/test_fun_hyper_dual.m 0 → 100644
+ 32
− 0
View file @ 2589646e
function test_fun_hyper_dual(f)
%
df = diff(f,'x');
ddf= diff(df,'x');
g = matlabFunction(f);
%
h1=0.000000001;
x = rand(4);
y = rand(4);
A = hyper_dual(x,ones(4)*h1, ones(4)*h1, zeros(4) );
ret = g(A,y);
if double(norm(ret.b / h1 - df(x,y),'fro')) < 10^-13 && ...
double(norm(ret.d / h1^2 - ddf(x,y),'fro')) < 10^-10 && ...
double(norm(ret.a - f(x,y),'fro')) < 10^-13
disp('hyper_dual ok')
elseif norm(ret.a - f(x,y),'fro') > 10^-13
disp('error in hyper_dual - evaluation')
disp(double(norm(ret.a - f(x,y),'fro')))
elseif norm(ret.b / h1 - df(x,y),'fro') > 10^-13
disp('error in hyper_dual - first derivative')
disp(double(norm(ret.b / h1 - df(x,y),'fro')))
else
disp('error in hyper_dual - second derivative')
disp(double(norm(ret.d / h1^2 - ddf(x,y),'fro')))
end
end
matlab/hyper_dual_numbers/test_hyper_dual.m 0 → 100644
+ 49
− 0
View file @ 2589646e
clear all
close all
whos
syms x y
%
disp('set 1')
f1 = symfun( x.^2 + y.^2-1 , [x y]);
test_fun_hyper_dual(f1)
%
f11 = symfun( 1 - x.^2 , [x y]);
test_fun_hyper_dual(f11)
%
f12 = symfun( x + y , [x y]);
test_fun_hyper_dual(f12)
%
f13 = symfun( -x .* y - y .* x , [x y]);
test_fun_hyper_dual(f13)
%
f14 = symfun( x .* y + y .* x , [x y]);
test_fun_hyper_dual(f14)
%
disp('set 2')
f2 = symfun( x ./ y , [x y]);
test_fun_hyper_dual(f2)
%
f21 = symfun( -x ./ y , [x y]);
test_fun_hyper_dual(f21)
%
f22 = symfun( -y ./ (x+1) , [x y]);
test_fun_hyper_dual(f22)
%
f23 = symfun( x ./ y - y ./ x , [x y]);
test_fun_hyper_dual(f23)
%
disp('set 3')
f3 = symfun( sin(x) .* cos(x) .* exp(x) , [x y]);
test_fun_hyper_dual(f3)
%
f31 = symfun( sinh(x) .*cosh(y.*x) .* sqrt(x.^y) , [x y]);
test_fun_hyper_dual(f31)
%
disp('set 4')
f41 = symfun( x.^3 , [x y]);
test_fun_hyper_dual(f41)
%
f42 = symfun( y .* 3.^x , [x y]);
test_fun_hyper_dual(f42)
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment