Add Matlab Programs

Matlab-Programs/PlotTest.m

+[t,y] = ode45(@cluster,[0:0.01:1],[1 2 3]);
+
+
+% figure(1)
+% plot(t,y(:,3));                 % plot of z(t) versus time
+% figure(2)
+% plot(t,y(:,1)); 
+% figure(3)
+% plot(y(:,1),y(:,3));            % plot of z versus x
+% figure(4)
+% plot3(y(:,1),y(:,2),y(:,3));    % 3D plot of trajectory
+% figure(5)
+% plot(y(:,1),y(:,2));            % plot of z versus x
+% figure(6)
+% plot(y(:,3),y(:,1));
+
+
+[x1,y1,z1] = meshgrid(-2:0.2:2,-2:0.2:2,-2:0.2:2);
+u = zeros(size(x1));
+v = zeros(size(y1));
+w = zeros(size(z1));
+
+t=0;
+for i = 1:numel(x1)
+    Yprime = cluster(t,[x1(i); y1(i); z1(i)]);
+    u(i) = Yprime(1);
+    v(i) = Yprime(2);
+    w(i) = Yprime(3);
+end
+
+
+quiver3(x1,y1,z1,u,v,w); figure(gcf)
\ No newline at end of file

Matlab-Programs/VanDePolOscillator.m

+
+
+
+mu = 100;
+
+
+F = @(t,y) [y(2); mu*(1-y(1)^2)*y(2)-y(1)];
+
+
+y0 = [0  1]';
+
+opts = odeset('stats','on')
+
+tspan= (0:1/36:1)*2*pi;
+
+tic
+[t,y] = ode15s(F,[0 460],y0,opts)    % stiff solver
+toc
+
+
+
+% plot(t,y(:,1), '.')
+% axis square
+% axis(1.2*[-1 1 -1 1])
+
+% plot phase-plane 
+plot(y(:,1),y(:,2),'.-')
\ No newline at end of file

Matlab-Programs/classifyMIN.m

+function [A, angle, type] = classifyMIN (mu_1,rho_1,a,b,t,set_mu_gamma,print_output)
+
+% returns
+%   A : Matrix of basis coefficients [a1,a2,a3]
+%
+%   type : 
+%   Type of minimizer 1 = (I) , 2 = (II) , 3 = (III) , 4 = (IV)
+%
+
+
+
+type = 0; % either 1,2,3,4 
+
+
+mu_h = @(b,t) mu_1.*(b./(t+(1-t).*b));  % harmonic mean
+mu_bar = @(b,t) mu_1.*((1-t)+t.*b);     % mu_bar
+
+if (set_mu_gamma == 'q1')
+  mu_gamma = @(b,t) mu_h(b,t);
+end
+if (set_mu_gamma == 'q2')
+  mu_gamma = @(b,t) mu_bar(b,t);
+end
+if (set_mu_gamma == 'm')
+  mu_gamma = @(b,t) 0.5*(mu_h(b,t) + mu_bar(b,t));
+end
+
+%  q1 q2 q3..
+q1 = mu_h(b,t);
+q2 = mu_bar(b,t); 
+q3 = mu_gamma(b,t);
+
+
+
+
+% values for q1,q2,q3 should be positiv
+% assert((q1 > 0 ) & (q2 > 0 ) & (q3 > 0), 'At least one of q1,q2 or q3 is not positive' )
+
+
+% Compute components of B_eff
+b1 = (mu_1*rho_1/4).*(b./(t+(1-t).*b)).*(1-t.*(1+a));
+b2 =  mu_1.*(rho_1/8).*(1-t.*(1+b.*a));
+
+
+
+% H = [q1 q3; q3 q2];          % right ??? 
+%check condition of H first
+% fprintf('condition number of Matrix H: %d \n', cond(H));
+
+
+
+
+epsilon = 1.e-18;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PARABOLIC CASE  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% if abs(det(A)) < epsilon * min(abs(det(A)),0)  
+if abs(q1*q2-q3^2) < epsilon     
+    
+    fprintf('determinant equal zero (parabolic case)')
+    fprintf('SHOULD NOT HAPPEN')
+    
+    % TODO 
+    if ( (b1/b2) - (q3/q1) < epsilon * min((b1/b2),(q3/q1)) )  
+        
+        %Minimizer not unique
+        type = 4;
+        % pick one arbitrary Minimzer as output...(TODO)
+         a1 = 1; %(TEST)
+         a2 = 2; %(TEST) 
+    
+    else
+        % TODO Hier weitere Fallunterscheidung nach b1b2 >0 ,< 0 nötig!
+        type = 4; % (TEST)
+        a1 = 0*b1; % (TEST)
+        a2 = b2; % (TEST)
+    end
+    
+end
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ELLIPTIC CASE  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+if (q1*q2-q3^2 > epsilon)
+%     fprintf('determinant greater than zero (elliptic case)');
+    
+%     a1_star =  (q2(b,t).*q1(b,t).*b1(a,b,t) - q3(b,t).*q2(b,t).*b2(a,b,t))./(q1(b,t).*q2(b,t)-q3(b,t).^2); 
+%     a2_star =  (q1(b,t).*q2(b,t).*b2(a,b,t) - q3(b,t).*q1(b,t).*b1(a,b,t))./(q1(b,t).*q2(b,t)-q3(b,t).^2); 
+    a1_star =  (q2.*q1.*b1 - q3.*q2.*b2)./(q1.*q2-q3.^2); 
+    a2_star =  (q1.*q2.*b2 - q3.*q1.*b1)./(q1.*q2-q3.^2);
+    
+    prod = a1_star*a2_star;
+    
+    
+    if(prod > epsilon) % (E1)  inside Lamba             %
+        % (a1_star,a2_star) is unique minimizer lies inside Lambda 
+        % therefore Minimizer not aligned with axes
+        
+%         fprintf('\n elliptic-case: (E1)');
+        a1 = a1_star;
+        a2 = a2_star;
+        type = 3;
+    end
+    % Make distinction between boundary & outside (prod < 0 ) 
+    if(abs(prod) < epsilon)  % (E2) on boundary of Lambda 
+%         fprintf('\n elliptic-case: (E2)');
+          
+        % Uniqueness of gloal minimizer if lies on boundary  if prod = 0 
+        % ----- % 
+        
+
+        % global minimizer lies on the boundary of Lambda depending on
+        % condition:
+        if (q2*b2^2 < q1*b1^2)  % global Minimizer given by (b1,0) 
+            a1 = b1;
+            a2 = 0*b1;
+            type = 1;  % Minimizer aligned with x1-axis
+        end
+        if  (q2*b2^2 > q1*b1^2)% global Minimizer given by (0,b2)    
+            a1 = 0*b1;
+            a2 = b2;
+            type = 2; % Minimizer aligned with x2-axis
+        end
+        if abs(q1*b1^2-q2*b2^2) < epsilon * min(q2*b2^2,q1*b1^2) % q1b1^2 = q2b2^2
+            % two Minimizers ..pick one
+            a1 = b1;
+            a2 = 0*b1;
+            type = 4;
+        end
+    end
+    if((prod) < -1*epsilon) %Outside of Lambda 
+        
+        if (q2*b2^2 < q1*b1^2)  % global Minimizer given by (b1,0) 
+            a1 = b1;
+            a2 = 0*b1;
+            type = 1;  % Minimizer aligned with x1-axis
+        end
+        if  (q2*b2^2 > q1*b1^2)% global Minimizer given by (0,b2)    
+            a1 = 0*b1;
+            a2 = b2;
+            type = 2; % Minimizer aligned with x2-axis
+        end
+    end
+    
+end
+    
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% HYPERBOLIC CASE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+if (q1*q2-q3^2 < -1*epsilon) 
+%     fprintf('determinant less than zero (hyperbolic case)');
+%     
+%     if (b1*b2 > 0) % (H1) 
+%         % can also be type (I) or (II) 
+%         %Minimizer either A1 = (0,b2) or A2 = (b1,0)
+%         % needs to be stable
+%         
+%         
+%         % check STABILITY 
+%         %check if A2 = (b1,0) is stable
+%         if (b1 > 0 && (q3*b1-q2*b2 > 0 )) 
+%             a1 = b1;
+%             a2 = 0*b2;
+%             type = 1;
+%         end
+%         if (b1 < 0 && (q3*b1-q2*b2 < 0 ) )
+%             a1 = b1;
+%             a2 = 0*b2;
+%             type = 1;
+%         end 
+%         %check if A1 = (0,b2) is stable 
+%          if (b2 > 0 && (q3*b2-q1*b1 > 0 ) )
+%             a1 = b1*0;
+%             a2 = b2;
+%             type = 2;
+%         end
+%         if (b2 < 0 && (q3*b2-q1*b1 < 0 )  )
+%             a1 = b1*0;
+%             a2 = b2;
+%             type = 2;
+%         end 
+% 
+%     end
+%     
+%     
+%     if ( abs(b1*b2) < epsilon )     %b1*b2 = 0 
+%         
+%         if (abs(b1) < epsilon )
+%             a1 = 0;
+%             a2 = b2;
+%         else % b2 = 0
+%             a1 = b1;
+%             a2 = 0; 
+%         end
+%         
+%     end
+% 
+%     if (b1*b2 < 0)  
+%   
+%         
+%         if (q2*b2^2 < q1*b1^2)  % global Minimizer given by (b1,0) 
+%             a1 = b1;
+%             a2 = 0*b2;
+%             type = 1;  % Minimizer aligned with x1-axis
+%         end
+%         if  (q2*b2^2 > q1*b1^2)% global Minimizer given by (0,b2)    
+%             a1 = 0*b1;
+%             a2 = b2;
+%             type = 2; % Minimizer aligned with x2-axis
+%         end
+%         if abs(q1*b1^2-q2*b2^2) < epsilon * min(q2*b2^2,q1*b1^2) % q1b1^2 = q2b2^2
+%             % two Minimizers ..pick one
+%             a1 = b1;
+%             a2 = 0*b1;
+%             type = 4;
+%         end
+%         
+%     end
+
+    if (q2*b2^2 < q1*b1^2)  % global Minimizer given by (b1,0) 
+        a1 = b1;
+        a2 = 0*b1;
+        type = 1;  % Minimizer aligned with x1-axis
+    end
+    if  (q2*b2^2 > q1*b1^2)% global Minimizer given by (0,b2)    
+        a1 = 0*b1;
+        a2 = b2;
+        type = 2; % Minimizer aligned with x2-axis
+    end
+    if abs(q1*b1^2-q2*b2^2) < epsilon * min(q2*b2^2,q1*b1^2) % q1b1^2 = q2b2^2
+        % two Minimizers ..pick one
+        a1 = b1;
+        a2 = 0*b1;
+        type = 4;
+    end
+
+
+    
+%     if ( abs(b1*b2) < epsilon)  % b1*b2 = 0 
+%     end
+    
+    
+    
+    % CAN NOT BE TYPE 3!!
+    
+end
+
+
+
+
+
+
+
+% Compute a3 from a1 % a2
+a3 = sqrt(2*a1*a2);
+
+
+
+
+% compute angle between [sqrt(a1) , sqrt(a2)] and e1:
+% angle = atan2(sqrt(a2),sqrt(a1));
+if (type == 3 )
+   angle = atan2(a2,a1);
+else
+   angle = 0;
+end
+% angle = atan2(norm(cross(a,b)), dot(a,b))
+
+
+%compute Kappa? 
+k = sqrt(abs(a1) + abs(a2));  % ? 
+
+
+
+
+% Coefficients of minimizer 
+
+if(print_output)
+   fprintf(' \n ') 
+   fprintf('=============================== OUTPUT ========================================= \n')
+   fprintf('Minimizer is of Type: %d  \n' , type);
+   fprintf('Coefficients a1,a2,a3 given by : %d, %d, %d   \n', a1, a2, a3);
+   fprintf('================================================================================ \n')
+end
+
+
+A = [a1, a2, a3];
+
+
+
+end
+

Matlab-Programs/cluster.m

+function f = cluster(t,y)
+  %BD
+  a=1;
+  b=1.2;
+  %equilibrium values
+  %c1 equilbrium value
+    f = zeros(size(y));
+    f(1) = -50*a*y(1)-b*y(1)+15*a*y(1)*y(2)+20*a*y(2)*y(3)+y(2)*b+9*a*y(2)^2+6*a*y(1)^2-60*a*y(2)-80*a*y(3)+24*a*y(2)*y(3)+16*a*y(3)^2;
+    f(2) = 10*a*y(1) - a*y(1)*y(2) -4*a*y(1)*y(3) -2*a*y(1)^2 -b*y(2) +3*a*y(2)^2 -10*a*y(2) +4*a*y(2)*y(3) +b*y(3);
+    f(3) = -2*a*y(1)*y(2) - 3*a*y(2)^2 -4*a*y(2)*y(3) +10*a*y(2)-b*y(3);
+end
+
+
+

Matlab-Programs/compute_F.m

+function F = compute_F(alpha,B,q1,q2,q3)
+
+% r = compute_r(alpha,B,q1,q2,q3);
+
+v = [cos(alpha);sin(alpha)];
+
+b1 = B(1,1);
+b2 = B(2,2);
+b3 = B(1,2);
+
+%compute Q(v_alpha x v_alpha)
+Q = q1.*v(1).^4 + q2.*v(2).^4 + 2.*q3.*v(1).^2.*v(2).^2;
+% 
+% TP = v*v';
+% L = stVenant(TP,mu,lambda);
+
+tmp1 = q1.*(v(1).^2+b1).^2 + q2.*(v(2).^2+b2).^2 + q3.*(sqrt(2)*v(1)*v(2)+b3).^2;
+tmp2 = q1.*b1.^2 + q2.*b2.^2+ q3.*b3.^2;
+L = 0.5*(tmp1-Q-tmp2) ;       %Polarization identity 
+
+
+r = L./Q;
+
+
+% F = r.^2.*Q - 2.*r.*trace(L'*B)
+F = (r.^2).*Q - 2.*r.*L
+
+
+end
+

Matlab-Programs/compute_r.m

+function r = compute_r(alpha,B,q1,q2,q3)
+
+v = [cos(alpha);sin(alpha)];
+
+b1 = B(1,1);
+b2 = B(2,2);
+b3 = B(1,2);
+
+
+% TP = v*v';
+% L = stVenant(TP,mu,lambda);
+
+
+
+%compute Q(v_alpha x v_alpha)
+Q = q1.*v(1).^4 + q2.*v(2).^4 + 2.*q3.*v(1).^2.*v(2).^2;
+
+tmp1 = q1.*(v(1).^2+b1).^2+q2.*(v(2).^2+b2).^2 + q3.*(sqrt(2)*v(1)*v(2)+b3).^2;
+tmp2 = q1.*b1.^2 + q2.*b2.^2+ q3.*b3.^2;
+L = 0.5.*(tmp1-Q-tmp2);         %Polarization identity 
+
+
+
+% r = trace(L'*B)./ Q;
+r = L./Q;
+
+
+end
+

Matlab-Programs/example1.m

+
+
+
+
+
+
+F = @(t,y) 2*y;
+
+t0 = 0;
+h = 1;
+tfinal = 3;
+y0 = 10;
+
+ode1(F,t0,h,tfinal,y0);
\ No newline at end of file

Matlab-Programs/exp.txt

+     x       exp(x)
+  0.00   1.00000000
+  0.10   1.10517092
+  0.20   1.22140276
+  0.30   1.34985881
+  0.40   1.49182470
+  0.50   1.64872127
+  0.60   1.82211880
+  0.70   2.01375271
+  0.80   2.22554093
+  0.90   2.45960311
+  1.00   2.71828183

Matlab-Programs/harmonicOscillator.m

+
+
+
+
+
+
+F = @(t,y) [y(2); -y(1)];
+
+
+y0 = [0  1]';
+
+
+
+tspan= (0:1/36:1)*2*pi;
+
+[t,y] = ode45(F,tspan,y0)
+
+
+% plot phase-plane 
+
+plot(y(:,1),y(:,2), 'o-')
+axis square
+axis(1.2*[-1 1 -1 1])
\ No newline at end of file

Matlab-Programs/isosurfaceTest.m

+[x,y,z] = meshgrid(1:20,1:20,1:20);
+data = sqrt(x.^2 + y.^2 + z.^2);
+p = patch(isosurface(x,y,z,data,20));
+isonormals(x,y,z,data,p)
+[r,g,b] = meshgrid(20:-1:1,1:20,1:20);
+isocolors(x,y,z,r/20,g/20,b/20,p)
+p.FaceColor = 'interp';
+p.EdgeColor = 'none';
+view(150,30) 
+daspect([1 1 1])
+camlight 
+lighting gouraud
\ No newline at end of file

Matlab-Programs/license.txt

+Copyright (c) 2009, Adam Auton
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are
+met:
+
+    * Redistributions of source code must retain the above copyright
+      notice, this list of conditions and the following disclaimer.
+    * Redistributions in binary form must reproduce the above copyright
+      notice, this list of conditions and the following disclaimer in
+      the documentation and/or other materials provided with the distribution
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
+LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+POSSIBILITY OF SUCH DAMAGE.