Sorted example scripts

MatlabGP/+ODE/+test/NS_1D_RHS.m

@@ -1,4 +1,4 @@
-function ydot = NS_1D_RHS(t,y,mu,K,gamma,dx)
+function ydot = NS_1D_RHS(~,y,mu,K,gamma,dx)
 
 if nargin<4
     mu = 1e-3;
@@ -13,10 +13,6 @@
 rhoU = y(:,2);
 E = y(:,3);
 
-% rho = utils.Filter(rho);
-% rhoU = utils.Filter(rhoU);
-% E = utils.Filter(E);
-
 rhoU(1) = 0;
 rhoU(end) = 0;
 
@@ -26,21 +22,15 @@
 rho(end) = rho(end-1);
 E(end) = E(end-1);
 
-
-
 U = rhoU./rho;
 
 dU = utils.Grad(U,dx);
-%dU = utils.Filter(dU);
 
 rhoU2 = rhoU.*U;
-%rhoU2 = utils.Filter(rhoU2);
 
 P = (gamma-1)*(E - 0.5*rhoU2);
-%P = utils.Filter(P);
 
 T = P./(R*rho);
-%T = utils.Filter(T);
 
 dT = utils.Grad(T,dx); 
 dT = utils.Filter(dT); 
@@ -49,18 +39,10 @@
 dF2 = utils.Grad(rhoU2 + P,dx);
 dF3 = utils.Grad((E+P).*U,dx);
 
-%dF1 = utils.Filter(dF1);
-%dF2 = utils.Filter(dF2);
-%dF3 = utils.Filter(dF3);
-
 dS2 = utils.Grad((4/3)*mu*dU,dx);
 dS3 = utils.Grad((4/3)*mu*U.*dU + K*dT,dx);
 
-%dS2 = utils.Filter(dS2);
-%dS3 = utils.Filter(dS3);
-
 ydot = [-dF1 dS2-dF2 dS3-dF3];
-%ydot = utils.Filter(ydot);
 
 ydot(1,:) = 0;
 ydot(end,:) = 0;

MatlabGP/+ODE/+test/NS_2D_RHS.m

@@ -0,0 +1,74 @@
+function ydot = NS_2D_RHS(~,y,mu,K,gamma,dx,dy)
+
+if nargin<3
+    mu = 1e-3;
+    K = 1e-3;
+    gamma = 1.4;
+    dx=1;
+    dy=1;
+end
+
+R = 8.314;
+
+rho = y(:,:,1);
+rhoU = y(:,:,2);
+rhoV = y(:,:,3);
+E = y(:,:,4);
+
+rhoU(1,:) = 0;
+rhoU(end,:) = 0;
+
+rhoV(1,:) = 0;
+rhoV(end,:) = 0;
+
+rho(1,:) = rho(2,:);
+E(1,:) = E(2,:);
+
+rho(end,:) = rho(end-1,:);
+E(end,:) = E(end-1,:);
+
+U = rhoU./rho;
+V = rhoV./rho;
+
+[dUy,dUx] = utils.Grad(U,dx,dy);
+[dVy,dVx] = utils.Grad(V,dx,dy);
+
+rhoU2 = rhoU.*U;
+rhoV2 = rhoV.*V;
+rhoUV = rhoU.*V;
+
+P = (gamma-1)*(E - 0.5*(rhoU2+rhoV2));
+
+T = P./(R*rho);
+
+[dTy,dTx] = utils.Grad(T,dx,dy); 
+
+dF1 = utils.Grad(rhoU,dx,dy);
+dF2 = utils.Grad(rhoU2 + P,dx);
+dF3 = utils.Grad(rhoUV,dx);
+dF4 = utils.Grad((E+P).*U,dx);
+
+[~,dG1] = utils.Grad(rhoV,dx,dy);
+[~,dG2] = utils.Grad(rhoUV,dx,dy);
+[~,dG3] = utils.Grad(rhoV2 + P,dx,dy);
+[~,dG4] = utils.Grad((E+P).*V,dx,dy);
+
+dS2x = utils.Grad((4/3)*mu*dUx,dx,dy);
+[~,dS2y] = utils.Grad((4/3)*mu*dUy,dx,dy);
+
+dS3x = utils.Grad((4/3)*mu*dVx,dx,dy);
+[~,dS3y] = utils.Grad((4/3)*mu*dVy,dx,dy);
+
+dS4x = utils.Grad((4/3)*mu*U.*dUx + K*dTx,dx,dy);
+[~,dS4y] = utils.Grad((4/3)*mu*V.*dVy + K*dTy,dx,dy);
+
+ydot(:,:,1) = -dF1 - dG1;
+ydot(:,:,2) = dS2x + dS2y -dF2 - dG2;
+ydot(:,:,3) = dS3x + dS3y -dF3 - dG3;
+ydot(:,:,4) = dS4x + dS4y -dF4 - dG4;
+
+
+ydot(1,:,:) = 0;
+ydot(end,:,:) = 0;
+
+end

MatlabGP/+ODE/rkf45.m

@@ -42,7 +42,7 @@
     Z = yn + 16/135*k1 + 6656/12825*k3 + 28561/56430*k4 - 9/50*k5 + ...
         2/55*k6;   
 
-    E = norm(Y-Z); % local error estimation    
+    E = norm(Y-Z,'fro'); % local error estimation    
     h = 0.9 * h * (TOL/E)^(1/5); 
 
     y{k} = Y;

MatlabGP/+optim/Adam.m

@@ -0,0 +1,94 @@
+classdef Adam
+
+    properties
+
+        iter = 0;
+
+        lb = [];
+        ub = [];
+
+        lr
+
+        wd
+
+        beta1
+        beta2
+
+        mt
+        vt
+    end
+
+    methods
+
+        function obj = Adam(x0,varargin)
+
+            input=inputParser;
+            input.KeepUnmatched=true;
+            input.PartialMatching=false;
+            input.addOptional('lb',[]);
+            input.addOptional('ub',[]);
+            input.addOptional('beta1',0.9);
+            input.addOptional('beta2',0.999);
+            input.addOptional('lr',0.1);
+            input.addOptional('wd',0);
+            input.parse(varargin{:})
+            in=input.Results;
+
+            obj.lb = in.lb;
+            obj.ub = in.ub;
+
+            obj.lr = in.lr;
+
+            obj.mt = 0*x0;
+            obj.vt = 0*x0;
+
+            obj.beta1 = in.beta1;
+            obj.beta2 = in.beta2;
+
+            obj.wd = in.wd;
+
+
+        end
+
+        function [obj,x] = step(obj,x,dF)
+
+            obj.iter = obj.iter + 1;
+
+            dF = dF + obj.wd.*x;
+
+            obj.mt = obj.beta1*obj.mt + (1 - obj.beta1)*dF;
+            obj.vt = obj.beta2*obj.vt + (1 - obj.beta2)*dF.^2;
+
+            mth = obj.mt./(1 - obj.beta1);
+            vth = obj.vt./(1 - obj.beta2);
+
+
+            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+            %%%%%%%%%%%%%% End Update Algorithm params %%%%%%%%%%%
+            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+            %update parameters
+            x = x - obj.lr*mth./(sqrt(vth) + eps );
+
+            %reflective upper bound
+            if ~isempty(obj.ub)
+                for jj = 1:length(x)
+                    if x(jj)>obj.ub(jj)
+                        x(jj)=obj.ub(jj) - 0.1*abs(obj.lr*mth(jj)./(sqrt(vth(jj)) + eps));
+                    end
+                end
+            end
+
+            %reflective lower bound
+            if ~isempty(obj.lb)
+                for jj = 1:length(x)
+                    if x(jj)<obj.lb(jj)
+                        x(jj)=obj.lb(jj) + 0.1*abs(obj.lr*mth(jj)./(sqrt(vth(jj)) + eps));
+                    end
+                end
+            end
+
+        end
+
+    end
+end

MatlabGP/+optim/SGD.m

@@ -0,0 +1,83 @@
+classdef SGD
+
+    %{
+    Stochastic Gradient Decent
+
+    input:
+    F - anonymous function to minimize (must return value and gradient)
+    x0 - initial guess point
+    
+    Optional Input:
+    lb - lower bound (reflective lower bound has been added)
+    ub - upper bound (reflective upper bound has been added)
+    lr - learning rate
+    iters - maximum number of iterations
+    tol - target tolerance on minimum
+
+    Output:
+    x - optimum location
+    Fx - value at optimum
+
+%}
+
+    properties
+
+        iter = 0;
+
+        lb = [];
+        ub = [];
+
+        lr
+
+        fv
+        xv
+    end
+
+    methods
+
+        function obj = SGD(x0,varargin)
+
+            input=inputParser;
+            input.KeepUnmatched=true;
+            input.PartialMatching=false;
+            input.addOptional('lb',[]);
+            input.addOptional('ub',[]);
+            input.addOptional('lr',0.1);
+            input.parse(varargin{:})
+            in=input.Results;
+
+            obj.lb = in.lb;
+            obj.ub = in.ub;
+
+            obj.lr = in.lr;
+        end
+
+        function [obj,x] = step(obj,x,dF)
+
+            obj.iter = obj.iter + 1;
+
+            %update parameters
+            x = x - obj.lr*dF;
+
+            %reflective upper bound
+            if ~isempty(obj.ub)
+                for jj = 1:length(x)
+                    if x(jj)>obj.ub(jj)
+                        x(jj)=obj.ub(jj) - 0.1*abs(obj.lr*dF(jj));
+                    end
+                end
+            end
+
+            %reflective lower bound
+            if ~isempty(obj.lb)
+                for jj = 1:length(x)
+                    if x(jj)<obj.lb(jj)
+                        x(jj)=obj.lb(jj) + 0.1*abs(obj.lr*dF(jj));
+                    end
+                end
+            end
+
+        end
+
+    end
+end