# Problem using function fsolve

46 views (last 30 days)
Vladyslav Toporin on 15 Jan 2020 at 17:42
Edited: Star Strider on 16 Jan 2020 at 21:14
Hello!
I'm trying to run a code with function fsolve.
T = 295.15;
p_a = 101.3;
x_0 = 0.07;
x_0_t = 0.05;
s_2 = 0.0225;
s_1 = 0.18;
d = 24.3e-3;
T = 295.15;
p_a = 101.3;
x_0 = 0.07;
x_0_t = 0.05;
s_2 = 0.0225;
s_1 = 0.18;
d = 24.3e-3;
rho_0=1.186
p_0=101.325;
c_0=343.2*sqrt(T/293);
Z_0 = rho_0*c_0;
Trennfrequenz=762;
M=real(Absorp_TransmissGreiner_RG2(1:100, 1:6));
F_=real(Absorp_TransmissGreiner_RG2_X(1:100));
lambda_0=c_0./F_;
rho=rho_0*p_a*T_0/(p_0*T);
k_0= 2*pi./lambda_0+sqrt(F_).*(1.94e-2/(c_0*d));
P_1=M(:,5);
P_3=M(:,2);
P_4=M(:,3);
P_6=M(:,6);
H_13=P_3./P_1;
H_46=P_6./P_4;
R_lp = (exp(k_0.*1i.*(x_0))-H_13.*exp(k_0.*1i.*(x_0+s_1)))./(exp(k_0.*(-1i).*x_0)-H_13.*exp(k_0.*(-1i).*(x_0+s_1)));
R_la = (exp(k_0.*1i.*(x_0_t))-H_46.*exp(k_0.*1i.*(x_0_t+s_1)))./(exp(k_0.*(-1i).*x_0_t)-H_13.*exp(k_0.*(-1i).*(x_0_t+s_1)));
Z_l_I = Z_0*(1+R_lp)./(1-R_lp);
Z_a_I = Z_0*(1+R_la)./(1-R_la);
M_II=real(Absorp_TransmissGreiner_RG2(1:100, 1:6));
P_1_II=M_II(:,5);
P_3_II=M_II(:,2);
P_4_II=M_II(:,3);
P_6_II=M_II(:,6);
H_13_II=P_3_II./P_1_II;
H_46_II=P_6_II./P_4_II;
R_lp_II = (exp(k_0.*1i.*(x_0))-H_13_II.*exp(k_0.*1i.*(x_0+s_1)))./(exp(k_0.*(-1i).*x_0)-H_13_II.*exp(k_0.*(-1i).*(x_0+s_1)));
R_la_II = (exp(k_0.*1i.*(x_0_t))-H_46_II.*exp(k_0.*1i.*(x_0_t+s_1)))./(exp(k_0.*(-1i).*x_0_t)-H_13_II.*exp(k_0.*(-1i).*(x_0_t+s_1)));
Z_l_II = Z_0*(1+R_lp_II)./(1-R_lp_II);
Z_a_II = Z_0*(1+R_la_II)./(1-R_la_II);
M_III=real(Absorp_TransmissGreiner_RG2(1:100, 1:6));
P_1_III=M_III(:,5);
P_3_III=M_III(:,2);
P_4_III=M_III(:,3);
P_6_III=M_III(:,6);
H_13_III=P_3_III./P_1_III;
H_46_III=P_6_III./P_4_III;
R_lp_III = (exp(k_0.*1i.*(x_0))-H_13_III.*exp(k_0.*1i.*(x_0+s_1)))./(exp(k_0.*(-1i).*x_0)-H_13_III.*exp(k_0.*(-1i).*(x_0+s_1)));
R_la_III = (exp(k_0.*1i.*(x_0_t))-H_46_III.*exp(k_0.*1i.*(x_0_t+s_1)))./(exp(k_0.*(-1i).*x_0_t)-H_13_III.*exp(k_0.*(-1i).*(x_0_t+s_1)));
Z_l_III = Z_0*(1+R_lp_III)./(1-R_lp_III);
Z_a_III = Z_0*(1+R_la_III)./(1-R_la_III);
z0 = [0;0;0];
F = @(z) [z(1)+(z(2)+T_0+Z_a_I)*z(3)/(z(2)+z(3)+T_0+Z_a_I)-Z_l_I;
z(1)+(z(2)+T_0+Z_a_II)*z(3)/(z(2)+z(3)+T_0+Z_a_II)-Z_l_II;
z(1)+(z(2)+T_0+Z_a_III)*z(3)/(z(2)+z(3)+T_0+Z_a_III)-Z_l_III];
z = fsolve(F,z0);
Z_1 = z(1,1);
Z_2 = z(1,2);
Z_3 = z(1,3);
...
And getting following errors (in Attachment as well):
Warning: Trust-region-dogleg algorithm of FSOLVE cannot handle non-square systems; using
> In fsolve (line 310)
Error using levenbergMarquardt (line 16)
Objective function is returning undefined values at initial point. fsolve cannot continue.
Error in fsolve (line 417)
levenbergMarquardt(funfcn,x,verbosity,options,defaultopt,f,JAC,caller, ...
z = fsolve(F,z0);

Star Strider on 15 Jan 2020 at 17:58
Use a different initial parameter estimate vector.
Try this:
z0 = rand(3,1);
Your code divides the parameters by each other, so the result will be NaN for the ‘z0’ your are currently using. They should be defined and finite for the ‘z0’ here.

Show 1 older comment
Star Strider on 15 Jan 2020 at 21:56
As always, my pleasure!
The fsolve function is a root finding function. I have no idea what files you are importing or what you are doing with them. However if you are doing curve fitting or some other related optimization problem, fsolve is not the correct oprimization function. There are much more appropriate functions for those purposes.
I can likely help you with that if I know what you want to do.
Vladyslav Toporin on 16 Jan 2020 at 14:47
This program is used in acoustics. Input data are arrays of acoustic pressure P_1, P_2, P_3 and P_4 for 4 microphones 26000x1 each and arrays of frequncy 26000x1 as well, from 1 to 26000.
Three pressure measurements with different conditions were completed, and three files are being loaded into the code.
Then from pressures coefficients for this system of nonlinear equations are calculated: Z_a_I, Z_a_II, Z_a_III, Z_l_I, Z_l_II, Z_l_III, which are arrays of the same size as pressures.
Variables z(1), z(2), z(3) are acoustic resistances, which are afterwards used for calculating transmission and absorption coefficients and bulding graphs dependent of frequency.
Here is the whole code:
T = 295.15;
p_a = 101.3;
x_0 = 0.07;
x_0_t = 0.05;
s_2 = 0.0225;
s_1 = 0.18;
d = 24.3e-3;
T = 295.15;
p_a = 101.3;
x_0 = 0.07;
x_0_t = 0.05;
s_2 = 0.0225;
s_1 = 0.18;
d = 24.3e-3;
rho_0=1.186
p_0=101.325;
c_0=343.2*sqrt(T/293);
Z_0 = rho_0*c_0;
Trennfrequenz=762;
M=real(Absorp_TransmissGreiner_RG2(1:100, 1:6));
F_=real(Absorp_TransmissGreiner_RG2_X(1:100));
lambda_0=c_0./F_;
rho=rho_0*p_a*T_0/(p_0*T);
k_0= 2*pi./lambda_0+sqrt(F_).*(1.94e-2/(c_0*d));
P_1=M(:,5);
P_3=M(:,2);
P_4=M(:,3);
P_6=M(:,6);
H_13=P_3./P_1;
H_46=P_6./P_4;
R_lp = (exp(k_0.*1i.*(x_0))-H_13.*exp(k_0.*1i.*(x_0+s_1)))./(exp(k_0.*(-1i).*x_0)-H_13.*exp(k_0.*(-1i).*(x_0+s_1)));
R_la = (exp(k_0.*1i.*(x_0_t))-H_46.*exp(k_0.*1i.*(x_0_t+s_1)))./(exp(k_0.*(-1i).*x_0_t)-H_13.*exp(k_0.*(-1i).*(x_0_t+s_1)));
Z_l_I = Z_0*(1+R_lp)./(1-R_lp);
Z_a_I = Z_0*(1+R_la)./(1-R_la);
M_II=real(Absorp_TransmissGreiner_RG2(1:100, 1:6));
P_1_II=M_II(:,5);
P_3_II=M_II(:,2);
P_4_II=M_II(:,3);
P_6_II=M_II(:,6);
H_13_II=P_3_II./P_1_II;
H_46_II=P_6_II./P_4_II;
R_lp_II = (exp(k_0.*1i.*(x_0))-H_13_II.*exp(k_0.*1i.*(x_0+s_1)))./(exp(k_0.*(-1i).*x_0)-H_13_II.*exp(k_0.*(-1i).*(x_0+s_1)));
R_la_II = (exp(k_0.*1i.*(x_0_t))-H_46_II.*exp(k_0.*1i.*(x_0_t+s_1)))./(exp(k_0.*(-1i).*x_0_t)-H_13_II.*exp(k_0.*(-1i).*(x_0_t+s_1)));
Z_l_II = Z_0*(1+R_lp_II)./(1-R_lp_II);
Z_a_II = Z_0*(1+R_la_II)./(1-R_la_II);
M_III=real(Absorp_TransmissGreiner_RG2(1:100, 1:6));
P_1_III=M_III(:,5);
P_3_III=M_III(:,2);
P_4_III=M_III(:,3);
P_6_III=M_III(:,6);
H_13_III=P_3_III./P_1_III;
H_46_III=P_6_III./P_4_III;
R_lp_III = (exp(k_0.*1i.*(x_0))-H_13_III.*exp(k_0.*1i.*(x_0+s_1)))./(exp(k_0.*(-1i).*x_0)-H_13_III.*exp(k_0.*(-1i).*(x_0+s_1)));
R_la_III = (exp(k_0.*1i.*(x_0_t))-H_46_III.*exp(k_0.*1i.*(x_0_t+s_1)))./(exp(k_0.*(-1i).*x_0_t)-H_13_III.*exp(k_0.*(-1i).*(x_0_t+s_1)));
Z_l_III = Z_0*(1+R_lp_III)./(1-R_lp_III);
Z_a_III = Z_0*(1+R_la_III)./(1-R_la_III);
z0 = [0;0;0];
F = @(z) [z(1)+(z(2)+T_0+Z_a_I)*z(3)/(z(2)+z(3)+T_0+Z_a_I)-Z_l_I;
z(1)+(z(2)+T_0+Z_a_II)*z(3)/(z(2)+z(3)+T_0+Z_a_II)-Z_l_II;
z(1)+(z(2)+T_0+Z_a_III)*z(3)/(z(2)+z(3)+T_0+Z_a_III)-Z_l_III];
z = fsolve(F,z0);
Z_1 = z(1,1);
Z_2 = z(1,2);
Z_3 = z(1,3);
A = 1+Z_1/Z_3;
B = Z_0*Z_1+Z_0*Z_2+Z_0*Z_1*Z_2/Z_3;
C = 1/Z_3*Z_0;
D = 1+Z_2/Z_3;
H = [A, B; C, D];
T_m = H*[1, 0; 1/1000*Z_0, 1];
Z_t = T_m(1,1)/T_m(2,1);
R = (Z_t/Z_0-1)/(Z_t/Z_0+1);
alpha = 1 - R^2;
figure (1)
subplot(2,1,1)
plot(F(F_<=Trennfrequenz),alpha(F_<=Trennfrequenz,1))
ylim([0,1])
xlim([0,7000])
xlabel('Frequenz [Hz]')
set(ti,'FontWeight','bold')
grid on
subplot(3,1,2)
plot(F(F_<=Trennfrequenz),abs(R(F_<=Trennfrequenz,1)))
ylim([0,1])
xlim([0,7000])
xlabel('Frequenz [Hz]')
set(ti,'FontWeight','bold')
grid on
Star Strider on 16 Jan 2020 at 14:56
It appears that you are doing parameter estimation.
In order to help, I need to know the data you have and the function you want to fit to them in order to estimate the parameters. (Acoustics is not an area of my expertise, so it may be necessary for you to explain some things that may not be readily apparent to me.)

Vladyslav Toporin on 16 Jan 2020 at 17:22
The data are .mat - files:
Thank you very much! You've already helped me a lot!
I don't know if you have a wish to read the whole explanation of my topic, but for any case i've described shortly measuring
principle and attached scientific publication which my method is based on.
The calculation method is used here is taken from this publication:
Unfortunately the button "Attachments'' isn't working.
I'm not sure if parameter estimation is done here.
For my measurements the acoustic tube (Kundt tube) is used.
Measures are completed with three different terminations: anechoic, echoic and open tube
My tube has 6 microphones(for low-frequent and 2 for high-frequent area).
In this code only low frequencies are handled.
Now i'm gonna describe the algorithm in detail:
T = 295.15; %Air temperature
p_a = 101.3; %Air pressure
x_0 = 0.07; %Distance between specimen and closest microphone (#2)
x_0_t = 0.05; %Distance between termination and closest microphone (#4)
s_2 = 0.0225; % s1 and s2: Distanve between microphones in each mic couple
s_1 = 0.18;
d = 24.3e-3; %Tube diameterT = 295.15;
rho_0=1.186 %Air density
p_0=101.325; %Air pressure inside tube
c_0=343.2*sqrt(T/293); %Air speed
Z_0 = rho_0*c_0; %Air impedance
Trennfrequenz=762; %Border frequence (border between low and high frequency region, in this program
%only low frequencies are measured
M=real(Absorp_TransmissGreiner_RG2(1:100, 1:6)); %pressure values of 6 microphone
F_=real(Absorp_TransmissGreiner_RG2_X(1:100)); %frequency array
lambda_0=c_0./F_; %sound wave length
rho=rho_0*p_a*T_0/(p_0*T); %calculated air density
k_0= 2*pi./lambda_0+sqrt(F_).*(1.94e-2/(c_0*d)); %wave number
P_1=M(:,5); %pressure only for first microphone
P_3=M(:,2); %pressure only for third microphone
P_4=M(:,3); %pressure only for fourth microphone
P_6=M(:,6); %pressure only for sixth microphone
H_13=P_3./P_1; %transfer function between 3rd and 1st microphone
H_46=P_6./P_4; %transfer function between 6th and 4th microphone
R_lp = (exp(k_0.*1i.*(x_0))-H_13.*exp(k_0.*1i.*(x_0+s_1)))./(exp(k_0.*(-1i).*x_0)-H_13.*exp(k_0.*(-1i).*(x_0+s_1))); %Reflection coefficient of the load
R_la = (exp(k_0.*1i.*(x_0_t))-H_46.*exp(k_0.*1i.*(x_0_t+s_1)))./(exp(k_0.*(-1i).*x_0_t)-H_13.*exp(k_0.*(-1i).*(x_0_t+s_1))); %Reflection coefficient of the termination
Z_l_I = Z_0*(1+R_lp)./(1-R_lp); %Impedance of everything that is behind the first tube section
Z_a_I = Z_0*(1+R_la)./(1-R_la); %Impedance of everything that is behind the termination
load ('E:\Masterprojekt\Messdaten\Absorp+TransmissGreiner_RG25_20mm_schallhart') %Now everything repeated for next measure data an so forth
M_II=real(Absorp_TransmissGreiner_RG2(1:100, 1:6));
P_1_II=M_II(:,5);
P_3_II=M_II(:,2);
P_4_II=M_II(:,3);
P_6_II=M_II(:,6);
H_13_II=P_3_II./P_1_II;
H_46_II=P_6_II./P_4_II;
R_lp_II = (exp(k_0.*1i.*(x_0))-H_13_II.*exp(k_0.*1i.*(x_0+s_1)))./(exp(k_0.*(-1i).*x_0)-H_13_II.*exp(k_0.*(-1i).*(x_0+s_1)));
R_la_II = (exp(k_0.*1i.*(x_0_t))-H_46_II.*exp(k_0.*1i.*(x_0_t+s_1)))./(exp(k_0.*(-1i).*x_0_t)-H_13_II.*exp(k_0.*(-1i).*(x_0_t+s_1)));
Z_l_II = Z_0*(1+R_lp_II)./(1-R_lp_II);
Z_a_II = Z_0*(1+R_la_II)./(1-R_la_II);
M_III=real(Absorp_TransmissGreiner_RG2(1:100, 1:6));
P_1_III=M_III(:,5);
P_3_III=M_III(:,2);
P_4_III=M_III(:,3);
P_6_III=M_III(:,6);
H_13_III=P_3_III./P_1_III;
H_46_III=P_6_III./P_4_III;
R_lp_III = (exp(k_0.*1i.*(x_0))-H_13_III.*exp(k_0.*1i.*(x_0+s_1)))./(exp(k_0.*(-1i).*x_0)-H_13_III.*exp(k_0.*(-1i).*(x_0+s_1)));
R_la_III = (exp(k_0.*1i.*(x_0_t))-H_46_III.*exp(k_0.*1i.*(x_0_t+s_1)))./(exp(k_0.*(-1i).*x_0_t)-H_13_III.*exp(k_0.*(-1i).*(x_0_t+s_1)));
Z_l_III = Z_0*(1+R_lp_III)./(1-R_lp_III);
Z_a_III = Z_0*(1+R_la_III)./(1-R_la_III);
z0 = [0;0;0];
F = @(z) [z(1)+(z(2)+T_0+Z_a_I)*z(3)/(z(2)+z(3)+T_0+Z_a_I)-Z_l_I;
z(1)+(z(2)+T_0+Z_a_II)*z(3)/(z(2)+z(3)+T_0+Z_a_II)-Z_l_II;
z(1)+(z(2)+T_0+Z_a_III)*z(3)/(z(2)+z(3)+T_0+Z_a_III)-Z_l_III];
z = fsolve(F,z0);
Z_1 = z(1,1);
Z_2 = z(1,2);
Z_3 = z(1,3);
A = 1+Z_1/Z_3; %A,B,C,D - coefficients for calculating matrixes for further finding
B = Z_0*Z_1+Z_0*Z_2+Z_0*Z_1*Z_2/Z_3; % absorption and reflection coefficients
C = 1/Z_3*Z_0;
D = 1+Z_2/Z_3;
H = [A, B; C, D];
T_m = H*[1, 0; 1/1000*Z_0, 1];
Z_t = T_m(1,1)/T_m(2,1);
R = (Z_t/Z_0-1)/(Z_t/Z_0+1);
alpha = 1 - R^2;
figure (1) %Plotting alpha and R
subplot(2,1,1)
plot(F(F_<=Trennfrequenz),alpha(F_<=Trennfrequenz,1))
ylim([0,1])
xlim([0,7000])
xlabel('Frequenz [Hz]')
set(ti,'FontWeight','bold')
grid on
subplot(3,1,2)
plot(F(F_<=Trennfrequenz),abs(R(F_<=Trennfrequenz,1)))
ylim([0,1])
xlim([0,7000])
xlabel('Frequenz [Hz]')
set(ti,'FontWeight','bold')
grid on

#### 1 Comment

Star Strider on 16 Jan 2020 at 17:28
As always, my pleasure!
This is far from my areas of expertise. However I may understand enough of it to suggest a reasonable approach to estimating the parameters (assuming I understand the parameters that need to be estimated).
EDIT — (16 Jan 2020 at 21:15)
It is theoretically possible to estimate the ‘A’, ‘B’, ‘C’, ‘D’ and ‘Z_0’ parameters in Eq 5. However, I need to know (1) if that is what you want to do, and (2) what the data files represent with respect to the diagram that accompanies Eq 5 (that I believe is Fig. 2).
I know essentially nothing about what you want to do. I will help as I can, however just now I do not know what the data files represent or what parameters you want to estimate, or what model you want to use to estimate them.