Skip to content
Snippets Groups Projects
Commit 46776e3d authored by Klaus Böhnlein's avatar Klaus Böhnlein
Browse files

Added Matlab Programs

parent 29e93dea
No related branches found
No related tags found
No related merge requests found
Showing
with 647 additions and 0 deletions
File added
clear all
clc
% --- Change PolynomialDisplayStyle ----
% sympref('PolynomialDisplayStyle','ascend');
% sympref('AbbreviateOutput',false);
syms f_plus(v1,v2,q1,q2,q3,b1,b2,b3)
assume( q1 > 0)
assume( q2 > 0)
assume( q3 > 0)
% assume(q3 == q1);
assume(q3 >= q1)
assume(q2 >= q3)
v = [v1; v2];
%should be sqrt(2) instead of 2!
f_plus(v1,v2,q1,q2,q3,b1,b2,b3) = q1*v1^4 + q2*v2^4+2*q3*v1^2*v2^2-2*(q1*b1*v1^2+ q2*b2*v2^2+sqrt(2)*q3*b3*v1*v2);
f_minus(v1,v2,q1,q2,q3,b1,b2,b3) = q1*v1^4 + q2*v2^4+2*q3*v1^2*v2^2+2*(q1*b1*v1^2+ q2*b2*v2^2+sqrt(2)*q3*b3*v1*v2);
% ---- Fix parameters
% f_plus = subs(f_plus,b3,0) % set b3
% f_plus = subs(f_plus,q3,q1)
% f_plus = subs(f_plus,q1,40);
% f_plus = subs(f_plus,q3,63.9);
% f_plus = subs(f_plus,q2,34.9);
% f_plus = subs(f_plus,b1,4);
% f_plus = subs(f_plus,b2,2.4);
% f_plus = subs(f_plus,b3,2.4);
% f_plus = subs(f_plus,q1,40);
% f_plus = subs(f_plus,q3,20);
% f_plus = subs(f_plus,q2,50);
% f_plus = subs(f_plus,b1,5);
% f_plus = subs(f_plus,b2,6);
% f_plus = subs(f_plus,b3,7);
f_plus = subs(f_plus,q1,40);
f_plus = subs(f_plus,q3,63.9);
f_plus = subs(f_plus,q2,70);
rndm1 = randi([1 20],1,1);
rndm2 = randi([1 20],1,1);
rndm3 = randi([1 20],1,1);
f_plus = subs(f_plus,b1,rndm1);
f_plus = subs(f_plus,b2,rndm2);
f_plus = subs(f_plus,b3,rndm3);
f_minus = subs(f_minus,q1,40);
f_minus = subs(f_minus,q3,63.9);
f_minus = subs(f_minus,q2,70);
f_minus = subs(f_minus,b1,rndm1);
f_minus = subs(f_minus,b2,rndm2);
f_minus = subs(f_minus,b3,rndm3);
% Compute Gradient
df_plusx = diff(f_plus,v1);
df_plusy = diff(f_plus,v2);
df_minusx = diff(f_minus,v1);
df_minusy = diff(f_minus,v2);
eq1 = df_plusx == 0;
eq2 = df_plusy == 0;
eqns = [eq1, eq2]
eq3 = df_minusx == 0;
eq4 = df_minusy == 0;
eqns_minus = [eq3, eq4]
S = solve(eqns,v1,v2,'MaxDegree' , 5, 'Real', true);
S_minus = solve(eqns_minus,v1,v2,'MaxDegree' , 5, 'Real', true);
% S = solve(eqns,v1,v2,'MaxDegree' , 5);
% S = solve(eqns,v1,v2,'MaxDegree' , 5, 'IgnoreAnalyticConstraints',true);
% S = solve(eqns,v1,v2,'MaxDegree' , 5, 'IgnoreAnalyticConstraints',true, 'Real', true);
% S = solve(eqns)
A = S.v1;
B = S.v2;
A_minus = S_minus.v1;
B_minus = S_minus.v2;
% A = simplify(A);
% B = simplify(B)
% S_minus = solve(eqns_minus,v1,v2);
% S_minus.v1
% S_minus.v2
%---------- TEST ---------------------
% r = subs(subs(df_plusx,v1,A(3)),v2,B(3));
fprintf('Testing equation grad(f) = 0 with stationary points')
double(subs(subs(df_plusx,v1,A(1)),v2,B(1)))
double(subs(subs(df_plusx,v1,A(2)),v2,B(2)))
double(subs(subs(df_plusx,v1,A(3)),v2,B(3)))
% ------------------------------------
fprintf('print stationary points of f_plus:')
double(A)
double(B)
fprintf('print stationary points of f_minus:')
double(A_minus)
double(B_minus)
% determine global Minimizer from stationary points:
fprintf('function values at stationary points:')
T = arrayfun(@(v1,v2) double(f_plus(v1,v2,q1,q2,q3,b1,b2,b3)),A,B,'UniformOutput', false)
T = cell2mat(T);
Min_plus = min(T, [], 'all')
% determine global Minimizer from stationary points:
fprintf('function values at stationary points:')
T_minus = arrayfun(@(v1,v2) double(f_minus(v1,v2,q1,q2,q3,b1,b2,b3)),A_minus,B_minus,'UniformOutput', false)
T_minus = cell2mat(T_minus);
Min_minus = min(T_minus, [], 'all')
globalMinimizerValue = min(Min_plus,Min_minus)
% Plot function
fsurf(@(x,y) f_plus(x,y,q1,q2,q3,b1,b2,b3))
hold on
plot3(A,B,T,'g*')
% view(90,0)
% view(2)
figure
fsurf(@(x,y) f_minus(x,y,q1,q2,q3,b1,b2,b3))
hold on
plot3(A_minus,B_minus,T,'g*')
%Write to txt-File
fileID = fopen('txt.txt','w');
fprintf(fileID,'%s' , latex(S.v1));
fclose(fileID);
File added
clear all
clc
% INPUT Parameters
mu_1 = 1;
rho_1 = 1;
theta= 0.3;
alpha = 2;
beta= 2;
fprintf('===========================================================================================================================');
fprintf(' Possible cases for global minimizers: (I) , (II) , (III) , (IV) ')
fprintf('============================================== INPUT PARAMETERS ===========================================================');
fprintf('mu_1: %d rho:1: %d theta: %d', mu_1,rho_1,theta );
fprintf('===========================================================================================================================');
% choose u_gamma = q3 to be either q1, q2 or something in between
% set_mu_gamma = 'q1';
set_mu_gamma = 'q2'; %(hyperbolic-case)
% set_mu_gamma = 'm'; % mean of q1,q2
% print_output = true;
print_output = false;
% Test for fixed value
% [A,angle,type] = classifyMIN(mu_1,rho_1,alpha,beta,theta,set_mu_gamma,print_output);
% PLOT
x = linspace(-20,20,45); %~alpha
y = linspace(1.5,40,45); %~beta
z = linspace(0.05,0.95,45); %~theta
[X1,Y1] = meshgrid(x,y);
[A_2D,angle_2D,V_2D] = arrayfun(@(a,b)classifyMIN(mu_1,rho_1,a,b,theta,set_mu_gamma,print_output),X1,Y1,'UniformOutput', false);
[X,Y,Z] = meshgrid(x,y,z);
[A,angle_3D,V_3D] = arrayfun(@(a,b,theta)classifyMIN(mu_1,rho_1,a,b,theta,set_mu_gamma,print_output),X,Y,Z,'UniformOutput', false);
color_3D = cell2mat(V_3D);
color_2D = cell2mat(V_2D);
angle_2D = cell2mat(angle_2D);
angle_3D = cell2mat(angle_3D);
A = cell2mat(A);
A_2D = cell2mat(A_2D);
X1 = reshape(X1,[],1);
Y1 = reshape(Y1,[],1);
X = reshape(X,[],1);
Y = reshape(Y,[],1);
Z = reshape(Z,[],1);
color_2D = reshape(color_2D,[],1);
color_3D = reshape(color_3D,[],1);
angle_2D = reshape(angle_2D,[],1);
angle_3D = reshape(angle_3D,[],1);
% Structure result depending on Type/Color
% V_2D = reshape(V_2D,[],1);
V_2DT1 = (color_2D == 1);
V_2DT2 = (color_2D == 2);
V_2DT3 = (color_2D == 3);
V_2DT1 = reshape(V_2DT1,[],1);
V_2DT2 = reshape(V_2DT2,[],1);
V_2DT3 = reshape(V_2DT3,[],1);
X1T1 = V_2DT1.*X1;
Y1T1 = V_2DT1.*Y1;
X1T2 = V_2DT2.*X1;
Y1T2 = V_2DT2.*Y1;
X1T3 = V_2DT3.*X1;
Y1T3 = V_2DT3.*Y1;
%%% 2D - Plot (fixed Theta)
cm_2D = redblue(3);
scatter(X1T1,Y1T1,[], 'MarkerFaceColor',[0 0 1 ], 'MarkerEdgeColor','black');
colormap(cm_2D)
hold on
scatter(X1T2,Y1T2,[], 'MarkerFaceColor',[1 0 0], 'MarkerEdgeColor','black');
hold on
% scatter(X1T3,Y1T3,[],'MarkerFaceColor',[0 0 1 ], 'MarkerEdgeColor','black');
scatter(X1T3,Y1T3,[], angle_2D, 'filled','MarkerEdgeColor','black');
colormap(cm_2D)
legend('Type 1', 'Type 2', 'Type 3')
title('Fixed Theta Plot')
xlabel('alpha')
ylabel('beta')
hold off
%%% 3D - Plot
V_3DT1 = (color_3D == 1);
V_3DT2 = (color_3D == 2);
V_3DT3 = (color_3D == 3);
V_3DT1 = reshape(V_3DT1,[],1);
V_3DT2 = reshape(V_3DT2,[],1);
V_3DT3 = reshape(V_3DT3,[],1);
XT1 = V_3DT1.*X;
YT1 = V_3DT1.*Y;
ZT1 = V_3DT1.*Z;
XT2 = V_3DT2.*X;
YT2 = V_3DT2.*Y;
ZT2 = V_3DT2.*Z;
XT3 = V_3DT3.*X;
YT3 = V_3DT3.*Y;
ZT3 = V_3DT3.*Z;
cm = redblue(90);
figure
%fixed Color
% scatter3(XT1,YT1,ZT1, [], 'MarkerFaceColor',[0.75 0 0],'MarkerEdgeColor', 'none');
%variing Color
scatter3(YT1,ZT1,XT1, [], 'MarkerFaceColor',[0 0 1 ], 'MarkerEdgeColor','black');
colormap(cm);
hold on
scatter3(YT2,ZT2,XT2, [], 'MarkerFaceColor',[1 0 0],'MarkerEdgeColor', 'none');
hold on
% scatter3(XT3,YT3,ZT3, [],'MarkerFaceColor',[0 0 1 ],'MarkerEdgeColor', 'none');
scatter3(YT3,ZT3,XT3, [], angle_3D, 'filled');
legend('Type 1', 'Type 2', 'Type 3')
title('Classification of Minimizers, theta:')
% xlabel('alpha')
% ylabel('beta')
% zlabel('theta')
xlabel('beta')
ylabel('theta')
zlabel('alpha')
File added
File added
mu_1 = 1;
rho_1 = 1;
theta = 0.5;
alpha = -0.5;
beta = 0.5;
% Compute components of B_eff
b1 = (mu_1*rho_1/4).*(beta./(theta+(1-theta).*beta)).*(1-theta.*(1+alpha));
b2 = mu_1.*(rho_1/8).*(1-theta.*(1+beta.*alpha));
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
% q1 q2 q3..
q1 = mu_h(beta,theta);
q2 = mu_bar(beta,theta);
fprintf('--------------------')
fprintf(' alpha:%d' , alpha)
fprintf(' beta:%d ' , beta)
fprintf(' theta:%d ', theta )
fprintf('-------------------- \n')
fprintf('q1*b1^2:')
q1*b1^2
fprintf('q2*b2^2:')
q2*b2^2
\ No newline at end of file
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];
%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')
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 (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
% 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
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
File added
File added
function c = redblue(m)
%REDBLUE Shades of red and blue color map
% REDBLUE(M), is an M-by-3 matrix that defines a colormap.
% The colors begin with bright blue, range through shades of
% blue to white, and then through shades of red to bright red.
% REDBLUE, by itself, is the same length as the current figure's
% colormap. If no figure exists, MATLAB creates one.
%
% For example, to reset the colormap of the current figure:
%
% colormap(redblue)
%
% See also HSV, GRAY, HOT, BONE, COPPER, PINK, FLAG,
% COLORMAP, RGBPLOT.
% Adam Auton, 9th October 2009
if nargin < 1, m = size(get(gcf,'colormap'),1); end
if (mod(m,2) == 0)
% From [0 0 1] to [1 1 1], then [1 1 1] to [1 0 0];
m1 = m*0.5;
r = (0:m1-1)'/max(m1-1,1);
g = r;
r = [r; ones(m1,1)];
g = [g; flipud(g)];
b = flipud(r);
else
% From [0 0 1] to [1 1 1] to [1 0 0];
m1 = floor(m*0.5);
r = (0:m1-1)'/max(m1,1);
g = r;
r = [r; ones(m1+1,1)];
g = [g; 1; flipud(g)];
b = flipud(r);
end
c = [r g b];
<?xml version="1.0" encoding="UTF-8"?>
<addonCore>
<label>Red Blue Colormap</label>
<version>1.0.0.0</version>
<type>zip</type>
<identifier>e5698820-4a80-11e4-9553-005056977bd0</identifier>
<createdBy name="Adam Auton"/>
<image>resources/screenshot.png</image>
</addonCore>
<?xml version="1.0" encoding="UTF-8"?>
<paths>
<path>.</path>
</paths>
<?xml version="1.0" encoding="UTF-8"?>
<addOn>
<identifier>e5698820-4a80-11e4-9553-005056977bd0</identifier>
<displayType>Function</displayType>
<translatedDisplayType>
<en_US>Function</en_US>
<ja_JP>関数</ja_JP>
<ko_KR>함수</ko_KR>
<zh_CN>函数</zh_CN>
</translatedDisplayType>
<name>Red Blue Colormap</name>
<author>Adam Auton</author>
<version>1.0.0.0</version>
<downloadUrl>https://www.mathworks.com/matlabcentral/mlc-downloads/downloads/submissions/25536/versions/1/download/zip?src=addons_ml_desktop_install&amp;profile_id=14169257&amp;license=40758619&amp;release_family=R2020b</downloadUrl>
<licenseUrl>https://addons.mathworks.com/registry/v1/e5698820-4a80-11e4-9553-005056977bd0/1.0.0.0/-/license</licenseUrl>
<previewImageUrl>https://www.mathworks.com/responsive_image/160/120/0/0/0/cache/matlabcentral/mlc-downloads/downloads/submissions/25536/versions/1/screenshot.png</previewImageUrl>
<releaseNotesUrl>https://www.mathworks.com/add-ons/e5698820-4a80-11e4-9553-005056977bd0/d1034754-ea44-08f8-b658-22b590dfae7e/releaseNotes</releaseNotesUrl>
<installationFolder>Functions</installationFolder>
</addOn>
Matlab-Programs/resources/previewImage.png

18.8 KiB

File added
Matlab-Programs/resources/screenshot.png

17 KiB

This diff is collapsed.
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment