Problem using function fsolve

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; 

load ('E:\Masterprojekt\Messdaten\Absorp+TransmissGreiner_RG25_20mm_offen')
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);  

load ('E:\Masterprojekt\Messdaten\Absorp+TransmissGreiner_RG25_20mm_schallhart')
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);             

load ('E:\Masterprojekt\Messdaten\Absorp+TransmissGreiner_RG25_20mm_Mit_Schaumstoff')
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

Levenberg-Marquardt algorithm instead.

> In fsolve (line 310)

In Absorptionsgrad_nach_Boonen_tiefe_Frequenz_korrigiert_1 (line 100)

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, ...

Error in Absorptionsgrad_nach_Boonen_tiefe_Frequenz_korrigiert_1 (line 100)

z = fsolve(F,z0);

Answer 1

[Star Strider]

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.

[Vladyslav Toporin]

Thank you for the help!

I tried your solution, but unfortunately the same error happened as before.

[Star Strider]

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]

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; 

load ('E:\Masterprojekt\Messdaten\Absorp+TransmissGreiner_RG25_20mm_offen')
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);  

load ('E:\Masterprojekt\Messdaten\Absorp+TransmissGreiner_RG25_20mm_schallhart')
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);     
       
load ('E:\Masterprojekt\Messdaten\Absorp+TransmissGreiner_RG25_20mm_Mit_Schaumstoff')
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]')
ylabel('Absorptionsgrad [-]')
ti=title('Absorptionsgrad');
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]')
ylabel('Reflexionsgrad [-]')
ti=title('Reflexionsgrad');
set(ti,'FontWeight','bold')
grid on

[Star Strider]

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.)

Answer 2

[Vladyslav Toporin]

The data are .mat - files:

https://www.dropbox.com/s/3gs7ubzerhjeiug/Absorp%2BTransmissGreiner_RG25_20mm_Mit_Schaumstoff.mat?dl=0

https://www.dropbox.com/s/9z7oj2hbm0mwtb5/Absorp%2BTransmissGreiner_RG25_20mm_offen.mat?dl=0

https://www.dropbox.com/s/rnmq5y8im89tk8x/Absorp%2BTransmissGreiner_RG25_20mm_schallhart.mat?dl=0

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:

http://pub.dega-akustik.de/ICA2019/data/articles/000567.pdf

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

load ('E:\Masterprojekt\Messdaten\Absorp+TransmissGreiner_RG25_20mm_offen') %loading measure data without termination
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);     
       
load ('E:\Masterprojekt\Messdaten\Absorp+TransmissGreiner_RG25_20mm_Mit_Schaumstoff')
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]')
ylabel('Absorptionsgrad [-]')
ti=title('Absorptionsgrad');
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]')
ylabel('Reflexionsgrad [-]')
ti=title('Reflexionsgrad');
set(ti,'FontWeight','bold')
grid on

[Star Strider]

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.