清空筛选条件
清空筛选条件

fminsearch and rcond with exponentials did not find expected results..i dont know if im doing something wrong or not,..please help

25 次查看（过去 30 天）
Odette
Odette 2024-6-18
评论： Odette 2024-6-18
Ran in:
clc;
c44 = 436+9;
e15 = -0.48;
e11 = 7.616-11;
u11 = 100;
d11 = 2.6;
p = 5700;
%c44e = 533.34;
e11e = 5.02e-11;
q = 1.602e-19;
uo = 1e20;
%i=sqrt(-1);
c44E = c44/e15;
E11E = e11 / e15;
d11E = d11 / (uo * u11);
h = 1;
m = uo * u11;
m1 = p * e11 * d11;
m2 = c44 * e11 * e15^2;
m3 = c44 * q * uo * u11;
m4 = uo * u11 * q;
a = d11 * m2;
omega_range = linspace(1, 10, 10); % Define your omega range
results = zeros(size(omega_range));
rcond_results = zeros(size(omega_range)); % Initialize array to store rcond results
for jj = 1:length(omega_range)
current_w = omega_range(jj);
% Define the objective function for fminsearch
objective_fn = @(k_vals) calculate_rcond(k_vals, current_w, h, c44, c44E, d11, e11, E11E, d11E, e15, m, m1, m2, m3, m4, a, p,e11e);
options=optimset('MaxIter',10000,'MaxFunEvals',1000);
% Initial guess for [k_real, k_imag]
initial_guess = [10, 10];
% Use fminsearch to find the values that minimize rcond
best_k = fminsearch(objective_fn, initial_guess,options);
% Store the best k value for the current omega
results(jj) = complex(best_k(1), best_k(2));
% Calculate rcond for the best_k found
rcond_A = calculate_rcond(best_k, current_w, h, c44, c44E, d11, e11, E11E, d11E, e15, m, m1, m2, m3, m4, a, p,e11e);
% Store the rcond value
rcond_results(jj) = rcond_A;
end
% Calculate phase velocity (omega / real(k))
c_real = omega_range ./ real(results);
% Plot omega vs. phase velocity
plot(omega_range, c_real, 'LineWidth', 2);
xlabel('\omega');
ylabel('c_real');
title('c_real vs. \omega');
% Display or plot the results as needed
disp('Results:');
Results:
disp(results);
1.0e+02 * Columns 1 through 9 7.0653 - 0.9996i 7.0713 - 8.1675i 7.0750 - 8.1692i 7.0750 - 8.1692i 7.0750 - 8.1692i 7.0754 - 0.9971i 7.0761 - 1.0031i 7.0784 - 1.0041i 7.0784 - 1.0041i Column 10 7.0766 + 1.3976i
% Check if there are roots (rcond close to zero)
tolerance = 1e-6; % Define a tolerance level for determining if rcond is "zero"
is_root = rcond_results < tolerance;
if any(is_root)
disp('Found roots for the determinant of the matrix A.');
disp('Corresponding omega and k values:');
for idx = find(is_root)
fprintf('Omega: %.2f, k: %.2f + %.2fi, rcond(A): %.2e\n', omega_range(idx), real(results(idx)), imag(results(idx)), rcond_results(idx));
end
else
disp('No roots found for the determinant of the matrix A within the specified tolerance.');
end
Found roots for the determinant of the matrix A.
Corresponding omega and k values:
Omega: 1.00, k: 706.53 + -99.96i, rcond(A): 0.00e+00 Omega: 2.00, k: 707.13 + -816.75i, rcond(A): 0.00e+00 Omega: 3.00, k: 707.50 + -816.92i, rcond(A): 0.00e+00 Omega: 4.00, k: 707.50 + -816.92i, rcond(A): 0.00e+00 Omega: 5.00, k: 707.50 + -816.92i, rcond(A): 0.00e+00 Omega: 6.00, k: 707.54 + -99.71i, rcond(A): 0.00e+00 Omega: 7.00, k: 707.61 + -100.31i, rcond(A): 0.00e+00 Omega: 8.00, k: 707.84 + -100.41i, rcond(A): 0.00e+00 Omega: 9.00, k: 707.84 + -100.41i, rcond(A): 0.00e+00 Omega: 10.00, k: 707.66 + 139.76i, rcond(A): 0.00e+00
function rcond_A = calculate_rcond(k_vals, current_w, h, c44, c44E, d11, e11, E11E, d11E, e15, m, m1, m2, m3, m4, a, p,e11e)
k_real = k_vals(1);
k_imag = k_vals(2);
% Calculate intermediate values b and c
b = m1 * current_w^2 / (k_real + k_imag*1i)^2 + current_w * m2 * 1i / (k_real + k_imag*1i)^2 - m3 / (k_real + k_imag*1i)^2;
c = (e11 * current_w * 1i - m4) * p * current_w^2 / (k_real + k_imag*1i)^4;
% Calculate s3 and s5 values for the current omega
s3_val = sqrt(1+(-b-sqrt(b^2-4*a * c)) /(2*a));
s5_val = sqrt(1+(-b+sqrt(b^2-4*a*c))/(2*a));
% Calculate p3, p5, q3, and q5 values for the current omega
p3_val = -e15*(s3_val^2-1)/(c44*(s3_val^2-1)+p*current_w^2/(k_real+k_imag*1i)^2);
p5_val = -e15*(s5_val^2-1)/(c44*(s5_val^2-1) + p*current_w^2/(k_real+k_imag*1i)^2);
q3_val = -m*(s3_val^2-1)/(d11 * (s3_val^2-1) + current_w*1i/ (k_real+k_imag*1i)^2);
q5_val = -m*(s5_val^2-1)/(d11 * (s5_val^2-1) + current_w*1i/ (k_real+k_imag*1i)^2);
% Define exp_kh and exp_minus_kh based on k_real and k_imag
exp_kh = exp((k_real+1i *k_imag)*h);
exp_minus_kh = exp(-(k_real+1i*k_imag)*h);
exp_s3=exp(s3_val*(k_real+1i*k_imag)*h);
exp_s4=exp(-s3_val*(k_real+1i*k_imag)*h);
exp_s5=exp(s5_val*(k_real+1i*k_imag)*h);
exp_s6=exp(-s5_val*(k_real+1i*k_imag)*h);
% Define matrix A
A = [exp_kh,exp_minus_kh,(c44E*p3_val+1)*s3_val*exp_s3,(c44E*p3_val+1)*(-s3_val)*exp_s4, ...
(c44E*p5_val+1)*s5_val*exp_s5,(c44E*p5_val+1)*(-s5_val)*exp_s6,0,0;
-E11E * exp_kh,E11E*exp_minus_kh, (p3_val-E11E)*s3_val * exp_s3, (p3_val - E11E) * (-s3_val) * exp_s4, ...
(p5_val-E11E) *s5_val * exp_s5, (p5_val - E11E) * (-s5_val) * exp_s6, 0, 0;
exp_kh,- exp_minus_kh, (1 - d11E * q3_val) * s3_val * exp_s3, (1 - d11E * q3_val) * (-s3_val) * exp_s4, ...
(1 - d11E * q5_val) * s5_val * exp_s5, (1 - d11E * q5_val) * (-s5_val) * exp_s6, 0, 0;
1, -1, (1 - d11E * q3_val) * s3_val, (1 - d11E * q3_val) * (-s3_val), ...
(1 - d11E * q5_val) * s5_val, (1 - d11E * q5_val) * (-s5_val), 0, 0;
1, -1, (c44E * p3_val + 1) * s3_val, (c44E * p3_val + 1) * (-s3_val), ...
(c44E * p5_val + 1) * s5_val, (c44E * p5_val + 1) * (-s5_val), (-c44E * s3_val) / e15, 0;
0, 0, p3_val, p3_val, p5_val, p5_val, -1, 0;
-E11E, E11E, (p3_val - E11E) * s3_val, (p3_val - E11E) * (-s3_val), ...
(p5_val - E11E) * s5_val, (p5_val - E11E) * (-s5_val), 0, e11e / e15;
1, 1, 1, 1, 1, 1, 0, -1];
% Compute the reciprocal condition number of the matrix A
rcond_A = rcond(A);
end

请先登录，再回答此问题。

回答（0 个）

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by