Erron on Symbolic calculation [Case 2] : "Empty sym : 0-by-1"

I need the help of people who have excellent talents in calculating using Matlab.

The equations to be implemented through Matlab coding is as follows.

Answer) β = 5.874395, A = 0.000060, b = 0.280599, φ = 5630.329851

The results to be obtained through this coding is written in "Answer) ~" above.

And, the data related to this problem(or eqution) is below.

I have tried with the following coding but kept getting the same results as the title. I would like to get helpful advice on this part.

I would appreciate it if you could give me a guide on the solution or error cause.

clear all
clc

th1 = [378 378 378 378; 0.4 0.4 0.4 0.4; 310 316 329 411]; % [V_i ; U_i ; T_i] 
th2 = [378 378 378 378; 0.8 0.8 0.8 0.8; 190 208 230 298];
th3 = [398 398 398 398; 0.4 0.4 0.4 0.4; 108 123 166 200];

syms B;
syms A;
syms phi;
syms b;

%------------- equation 1 -------------%
a1_th1 = 0;					
a1_th2 = 0;
a1_th3 = 0;

i = 1;					% "Stress Lv. 1" in the second term of ∂Λ/∂β(equation 1)
 for j = 1:(numel(th1(i, :)))
  N = numel(th1(i,j));			
                                                	
  a1_th1 = N*(1-(th1(i+2,j)/A*exp(-(phi/th1(i,j)+b/th1(i+1,j))))^B)*(log(th1(i+2,j)/A)-(phi/th1(i,j)+b/th1(i+1,j))) + a1_th1 ;
	
 end

i = 1;					% "Stress Lv. 2" in the second term of ∂Λ/∂β(equation 1)
 for j = 1:(numel(th2(i, :)))
  N = numel(th2(i,j));			
                                                	
  a1_th2 = N*(1-(th2(i+2,j)/A*exp(-(phi/th2(i,j)+b/th2(i+1,j))))^B)*(log(th2(i+2,j)/A)-(phi/th2(i,j)+b/th2(i+1,j))) + a1_th2 ;
	
 end

i = 1;					% "Stress Lv. 3" in the second term of ∂Λ/∂β(equation 1)
 for j = 1:(numel(th3(i, :)))
  N = numel(th3(i,j));			
                                                	
  a1_th3 = N*(1-(th3(i+2,j)/A*exp(-(phi/th3(i,j)+b/th3(i+1,j))))^B)*(log(th3(i+2,j)/A)-(phi/th3(i,j)+b/th3(i+1,j))) + a1_th3 ;
	
 end

a1 = a1_th1 + a1_th2 + a1_th3;
N = numel(th1(1,:)) + numel(th2(1,:)) + numel(th3(1,:));		% ∂Λ/∂β first term

one = N/B + a1;			
% vpa(one,3);


%------------- equation 2 -------------%
b1 = 0;
b2 = 0;
b3 = 0;

i = 1;					% "Stress Lv. 1" in the second term of ∂Λ/∂A(equation 2)
 for j = 1:(numel(th1(i, :)))
  N = numel(th1(i,j));
  b1 = N*(th1(i+2,j)/A*exp(-(phi/th1(i,j)+b/th1(i+1,j))))^B+b1;	
 end

i = 1;					% "Stress Lv. 2" in the second term of ∂Λ/∂A(equation 2)
 for j = 1:(numel(th2(i, :)))
  N = numel(th2(i,j));
  b2 = N*(th2(i+2,j)/A*exp(-(phi/th2(i,j)+b/th2(i+1,j))))^B+b2;	
 end

i = 1;					% "Stress Lv. 3" in the second term of ∂Λ/∂A(equation 2)
 for j = 1:(numel(th3(i, :)))
  N = numel(th3(i,j));
  b3 = N*(th3(i+2,j)/A*exp(-(phi/th3(i,j)+b/th3(i+1,j))))^B+b3;	
 end

N = numel(th1(1,:)) + numel(th2(1,:)) + numel(th3(1,:));		% ∂Λ/∂A first term

two = -B/A*(N + (b1 + b2 + b3));		  	
% vpa(two,3);

%------------- equation 3 -------------%
c1 = 0;
c2 = 0;
c3 = 0;

i = 1;					% "Stress Lv. 1" in the second term of ∂Λ/∂φ(equation 3)
 for j = 1:(numel(th1(i, :)))
  N = numel(th1(i,j));			
  c1 = N/th1(i,j)*(th1(i+2,j)/A*exp(-(phi/th1(i,j)+b/th1(i+1,j))))^B+c1;
 end

i = 1;					% "Stress Lv. 2" in the second term of ∂Λ/∂φ(equation 3)
 for j = 1:(numel(th2(i, :)))
  N = numel(th2(i,j));			
  c2 = N/th2(i,j)*(th2(i+2,j)/A*exp(-(phi/th2(i,j)+b/th2(i+1,j))))^B+c2;
 end

i = 1;					% "Stress Lv. 3" in the second term of ∂Λ/∂φ(equation 3)
 for j = 1:(numel(th3(i, :)))
  N = numel(th3(i,j));			
  c3 = N/th3(i,j)*(th3(i+2,j)/A*exp(-(phi/th3(i,j)+b/th3(i+1,j))))^B+c3;
 end

N_c1 = 0;				
N_c2 = 0;
N_c3 = 0;

i = 1;
 for j = 1:(numel(th1(i, :)))
  N_c1 = numel(th1(i,j))/th1(i,j)+N_c1;		% "Stress Lv. 1" in the first term of ∂Λ/∂φ(equation 3)
 end					

i = 1;
 for j = 1:(numel(th2(i, :)))
  N_c2 = numel(th2(i,j))/th2(i,j)+N_c2;		% "Stress Lv. 2" in the first term of ∂Λ/∂φ(equation 3)
 end	

i = 1;
 for j = 1:(numel(th3(i, :)))
  N_c3 = numel(th3(i,j))/th3(i,j)+N_c3;		% "Stress Lv. 3" in the first term of ∂Λ/∂φ(equation 3)
 end	

N_c = N_c1 + N_c2 + N_c3;				  % ∂Λ/∂φ first term

three = -B*(N_c  - (c1 + c2 + c3));		  	
% vpa(three,3);

%------------- equation 3 -------------%
d1 = 0;
d2 = 0;
d3 = 0;

i = 1;					% "Stress Lv. 1" in the second term of ∂Λ/∂b(equation 4)
 for j = 1:(numel(th1(i, :)))
  N = numel(th1(i,j));			
  d1 = N/th1(i+1,j)*(th1(i+2,j)/A*exp(-(phi/th1(i,j)+b/th1(i+1,j))))^B+d1;
 end

i = 1;					% "Stress Lv. 2" in the second term of ∂Λ/∂b(equation 4)"
 for j = 1:(numel(th2(i, :)))
  N = numel(th2(i,j));			
  d2 = N/th2(i+1,j)*(th2(i+2,j)/A*exp(-(phi/th2(i,j)+b/th2(i+1,j))))^B+d2;
 end

i = 1;					% "Stress Lv. 3" in the second term of ∂Λ/∂b(equation 4)
 for j = 1:(numel(th3(i, :)))
  N = numel(th3(i,j));			
  d3 = N/th3(i+1,j)*(th3(i+2,j)/A*exp(-(phi/th3(i,j)+b/th3(i+1,j))))^B+d3;
 end

N_d1 = 0;				
N_d2 = 0;
N_d3 = 0;

i = 1;
 for j = 1:(numel(th1(i, :)))
  N_d1 = numel(th1(i,j))/th1(i+1,j)+N_d1;		% "Stress Lv. 1" in the first term of ∂Λ/∂b(equation 4)
 end					

i = 1;
 for j = 1:(numel(th2(i, :)))
  N_d2 = numel(th2(i,j))/th2(i+1,j)+N_d2;		% "Stress Lv. 2" in the first term of ∂Λ/∂b(equation 4)
 end	

i = 1;
 for j = 1:(numel(th3(i, :)))
  N_d3 = numel(th3(i,j))/th3(i+1,j)+N_d3;		% "Stress Lv. 3" in the first term of ∂Λ/∂b(equation 4)
 end	

N_d = N_d1 + N_d2 + N_d3;				  % ∂Λ/∂b first term

four = -B*(N_d - (d1 + d2 + d3));		  
% vpa(four,3);

[B, A, phi, b] = solve(one,two,three,four,'Real',true)	 

% eqns = [one == 0, two == 0, three == 0, four == 0];
% vars = [B, A, phi, b];
% [solB, solA, solphi, solb] = solve(eqns, vars, 'Real', true)

[Star Strider]

If you display ‘one’, ‘two’, ‘three’, and ‘four’ you see the problem.

There may not be an analytic (symbolic) solution —

clear all
clc

th1 = [378 378 378 378; 0.4 0.4 0.4 0.4; 310 316 329 411]; % [V_i ; U_i ; T_i] 
th2 = [378 378 378 378; 0.8 0.8 0.8 0.8; 190 208 230 298];
th3 = [398 398 398 398; 0.4 0.4 0.4 0.4; 108 123 166 200];

syms B;
syms A;
syms phi;
syms b;

%------------- equation 1 -------------%
a1_th1 = 0;					
a1_th2 = 0;
a1_th3 = 0;

i = 1;					% "Stress Lv. 1" in the second term of ∂Λ/∂β(equation 1)
 for j = 1:(numel(th1(i, :)))
  N = numel(th1(i,j));			
                                                	
  a1_th1 = N*(1-(th1(i+2,j)/A*exp(-(phi/th1(i,j)+b/th1(i+1,j))))^B)*(log(th1(i+2,j)/A)-(phi/th1(i,j)+b/th1(i+1,j))) + a1_th1 ;
	
 end

i = 1;					% "Stress Lv. 2" in the second term of ∂Λ/∂β(equation 1)
 for j = 1:(numel(th2(i, :)))
  N = numel(th2(i,j));			
                                                	
  a1_th2 = N*(1-(th2(i+2,j)/A*exp(-(phi/th2(i,j)+b/th2(i+1,j))))^B)*(log(th2(i+2,j)/A)-(phi/th2(i,j)+b/th2(i+1,j))) + a1_th2 ;
	
 end

i = 1;					% "Stress Lv. 3" in the second term of ∂Λ/∂β(equation 1)
 for j = 1:(numel(th3(i, :)))
  N = numel(th3(i,j));			
                                                	
  a1_th3 = N*(1-(th3(i+2,j)/A*exp(-(phi/th3(i,j)+b/th3(i+1,j))))^B)*(log(th3(i+2,j)/A)-(phi/th3(i,j)+b/th3(i+1,j))) + a1_th3 ;
	
 end

a1 = a1_th1 + a1_th2 + a1_th3;
N = numel(th1(1,:)) + numel(th2(1,:)) + numel(th3(1,:));		% ∂Λ/∂β first term

one = N/B + a1;			
% vpa(one,3);


%------------- equation 2 -------------%
b1 = 0;
b2 = 0;
b3 = 0;

i = 1;					% "Stress Lv. 1" in the second term of ∂Λ/∂A(equation 2)
 for j = 1:(numel(th1(i, :)))
  N = numel(th1(i,j));
  b1 = N*(th1(i+2,j)/A*exp(-(phi/th1(i,j)+b/th1(i+1,j))))^B+b1;	
 end

i = 1;					% "Stress Lv. 2" in the second term of ∂Λ/∂A(equation 2)
 for j = 1:(numel(th2(i, :)))
  N = numel(th2(i,j));
  b2 = N*(th2(i+2,j)/A*exp(-(phi/th2(i,j)+b/th2(i+1,j))))^B+b2;	
 end

i = 1;					% "Stress Lv. 3" in the second term of ∂Λ/∂A(equation 2)
 for j = 1:(numel(th3(i, :)))
  N = numel(th3(i,j));
  b3 = N*(th3(i+2,j)/A*exp(-(phi/th3(i,j)+b/th3(i+1,j))))^B+b3;	
 end

N = numel(th1(1,:)) + numel(th2(1,:)) + numel(th3(1,:));		% ∂Λ/∂A first term

two = -B/A*(N + (b1 + b2 + b3));		  	
% vpa(two,3);

%------------- equation 3 -------------%
c1 = 0;
c2 = 0;
c3 = 0;

i = 1;					% "Stress Lv. 1" in the second term of ∂Λ/∂φ(equation 3)
 for j = 1:(numel(th1(i, :)))
  N = numel(th1(i,j));			
  c1 = N/th1(i,j)*(th1(i+2,j)/A*exp(-(phi/th1(i,j)+b/th1(i+1,j))))^B+c1;
 end

i = 1;					% "Stress Lv. 2" in the second term of ∂Λ/∂φ(equation 3)
 for j = 1:(numel(th2(i, :)))
  N = numel(th2(i,j));			
  c2 = N/th2(i,j)*(th2(i+2,j)/A*exp(-(phi/th2(i,j)+b/th2(i+1,j))))^B+c2;
 end

i = 1;					% "Stress Lv. 3" in the second term of ∂Λ/∂φ(equation 3)
 for j = 1:(numel(th3(i, :)))
  N = numel(th3(i,j));			
  c3 = N/th3(i,j)*(th3(i+2,j)/A*exp(-(phi/th3(i,j)+b/th3(i+1,j))))^B+c3;
 end

N_c1 = 0;				
N_c2 = 0;
N_c3 = 0;

i = 1;
 for j = 1:(numel(th1(i, :)))
  N_c1 = numel(th1(i,j))/th1(i,j)+N_c1;		% "Stress Lv. 1" in the first term of ∂Λ/∂φ(equation 3)
 end					

i = 1;
 for j = 1:(numel(th2(i, :)))
  N_c2 = numel(th2(i,j))/th2(i,j)+N_c2;		% "Stress Lv. 2" in the first term of ∂Λ/∂φ(equation 3)
 end	

i = 1;
 for j = 1:(numel(th3(i, :)))
  N_c3 = numel(th3(i,j))/th3(i,j)+N_c3;		% "Stress Lv. 3" in the first term of ∂Λ/∂φ(equation 3)
 end	

N_c = N_c1 + N_c2 + N_c3;				  % ∂Λ/∂φ first term

three = -B*(N_c  - (c1 + c2 + c3));		  	
% vpa(three,3);

%------------- equation 3 -------------%
d1 = 0;
d2 = 0;
d3 = 0;

i = 1;					% "Stress Lv. 1" in the second term of ∂Λ/∂b(equation 4)
 for j = 1:(numel(th1(i, :)))
  N = numel(th1(i,j));			
  d1 = N/th1(i+1,j)*(th1(i+2,j)/A*exp(-(phi/th1(i,j)+b/th1(i+1,j))))^B+d1;
 end

i = 1;					% "Stress Lv. 2" in the second term of ∂Λ/∂b(equation 4)"
 for j = 1:(numel(th2(i, :)))
  N = numel(th2(i,j));			
  d2 = N/th2(i+1,j)*(th2(i+2,j)/A*exp(-(phi/th2(i,j)+b/th2(i+1,j))))^B+d2;
 end

i = 1;					% "Stress Lv. 3" in the second term of ∂Λ/∂b(equation 4)
 for j = 1:(numel(th3(i, :)))
  N = numel(th3(i,j));			
  d3 = N/th3(i+1,j)*(th3(i+2,j)/A*exp(-(phi/th3(i,j)+b/th3(i+1,j))))^B+d3;
 end


 
N_d1 = 0;				
N_d2 = 0;
N_d3 = 0;

i = 1;
 for j = 1:(numel(th1(i, :)))
  N_d1 = numel(th1(i,j))/th1(i+1,j)+N_d1;		% "Stress Lv. 1" in the first term of ∂Λ/∂b(equation 4)
 end					

i = 1;
 for j = 1:(numel(th2(i, :)))
  N_d2 = numel(th2(i,j))/th2(i+1,j)+N_d2;		% "Stress Lv. 2" in the first term of ∂Λ/∂b(equation 4)
 end	

i = 1;
 for j = 1:(numel(th3(i, :)))
  N_d3 = numel(th3(i,j))/th3(i+1,j)+N_d3;		% "Stress Lv. 3" in the first term of ∂Λ/∂b(equation 4)
 end	

N_d = N_d1 + N_d2 + N_d3;				  % ∂Λ/∂b first term

four = -B*(N_d - (d1 + d2 + d3));		  
% vpa(four,3);

one

one = 

two

two = 

●

three

three = 

four

four = 

[B, A, phi, b] = solve(one,two,three,four,'Real',true)	 

Warning: Unable to find explicit solution. For options, see help.

 
B =
 
Empty sym: 0-by-1
 
 
A =
 
Empty sym: 0-by-1
 
 
phi =
 
Empty sym: 0-by-1
 
 
b =
 
Empty sym: 0-by-1
 

% eqns = [one == 0, two == 0, three == 0, four == 0];
% vars = [B, A, phi, b];
% [solB, solA, solphi, solb] = solve(eqns, vars, 'Real', true)

I am not even sure that a numeric solution would be possible.

.

[Walter Roberson]

Let's reformulate it as a minimization:

clear all
clc

th1 = [378 378 378 378; 0.4 0.4 0.4 0.4; 310 316 329 411]; % [V_i ; U_i ; T_i] 
th2 = [378 378 378 378; 0.8 0.8 0.8 0.8; 190 208 230 298];
th3 = [398 398 398 398; 0.4 0.4 0.4 0.4; 108 123 166 200];

syms B;
syms A;
syms phi;
syms b;

%------------- equation 1 -------------%
a1_th1 = 0;					
a1_th2 = 0;
a1_th3 = 0;

i = 1;					% "Stress Lv. 1" in the second term of ∂Λ/∂β(equation 1)
 for j = 1:(numel(th1(i, :)))
  N = numel(th1(i,j));			
                                                	
  a1_th1 = N*(1-(th1(i+2,j)/A*exp(-(phi/th1(i,j)+b/th1(i+1,j))))^B)*(log(th1(i+2,j)/A)-(phi/th1(i,j)+b/th1(i+1,j))) + a1_th1 ;
	
 end

i = 1;					% "Stress Lv. 2" in the second term of ∂Λ/∂β(equation 1)
 for j = 1:(numel(th2(i, :)))
  N = numel(th2(i,j));			
                                                	
  a1_th2 = N*(1-(th2(i+2,j)/A*exp(-(phi/th2(i,j)+b/th2(i+1,j))))^B)*(log(th2(i+2,j)/A)-(phi/th2(i,j)+b/th2(i+1,j))) + a1_th2 ;
	
 end

i = 1;					% "Stress Lv. 3" in the second term of ∂Λ/∂β(equation 1)
 for j = 1:(numel(th3(i, :)))
  N = numel(th3(i,j));			
                                                	
  a1_th3 = N*(1-(th3(i+2,j)/A*exp(-(phi/th3(i,j)+b/th3(i+1,j))))^B)*(log(th3(i+2,j)/A)-(phi/th3(i,j)+b/th3(i+1,j))) + a1_th3 ;
	
 end

a1 = a1_th1 + a1_th2 + a1_th3;
N = numel(th1(1,:)) + numel(th2(1,:)) + numel(th3(1,:));		% ∂Λ/∂β first term

one = N/B + a1;			
% vpa(one,3);


%------------- equation 2 -------------%
b1 = 0;
b2 = 0;
b3 = 0;

i = 1;					% "Stress Lv. 1" in the second term of ∂Λ/∂A(equation 2)
 for j = 1:(numel(th1(i, :)))
  N = numel(th1(i,j));
  b1 = N*(th1(i+2,j)/A*exp(-(phi/th1(i,j)+b/th1(i+1,j))))^B+b1;	
 end

i = 1;					% "Stress Lv. 2" in the second term of ∂Λ/∂A(equation 2)
 for j = 1:(numel(th2(i, :)))
  N = numel(th2(i,j));
  b2 = N*(th2(i+2,j)/A*exp(-(phi/th2(i,j)+b/th2(i+1,j))))^B+b2;	
 end

i = 1;					% "Stress Lv. 3" in the second term of ∂Λ/∂A(equation 2)
 for j = 1:(numel(th3(i, :)))
  N = numel(th3(i,j));
  b3 = N*(th3(i+2,j)/A*exp(-(phi/th3(i,j)+b/th3(i+1,j))))^B+b3;	
 end

N = numel(th1(1,:)) + numel(th2(1,:)) + numel(th3(1,:));		% ∂Λ/∂A first term

two = -B/A*(N + (b1 + b2 + b3));		  	
% vpa(two,3);

%------------- equation 3 -------------%
c1 = 0;
c2 = 0;
c3 = 0;

i = 1;					% "Stress Lv. 1" in the second term of ∂Λ/∂φ(equation 3)
 for j = 1:(numel(th1(i, :)))
  N = numel(th1(i,j));			
  c1 = N/th1(i,j)*(th1(i+2,j)/A*exp(-(phi/th1(i,j)+b/th1(i+1,j))))^B+c1;
 end

i = 1;					% "Stress Lv. 2" in the second term of ∂Λ/∂φ(equation 3)
 for j = 1:(numel(th2(i, :)))
  N = numel(th2(i,j));			
  c2 = N/th2(i,j)*(th2(i+2,j)/A*exp(-(phi/th2(i,j)+b/th2(i+1,j))))^B+c2;
 end

i = 1;					% "Stress Lv. 3" in the second term of ∂Λ/∂φ(equation 3)
 for j = 1:(numel(th3(i, :)))
  N = numel(th3(i,j));			
  c3 = N/th3(i,j)*(th3(i+2,j)/A*exp(-(phi/th3(i,j)+b/th3(i+1,j))))^B+c3;
 end

N_c1 = 0;				
N_c2 = 0;
N_c3 = 0;

i = 1;
 for j = 1:(numel(th1(i, :)))
  N_c1 = numel(th1(i,j))/th1(i,j)+N_c1;		% "Stress Lv. 1" in the first term of ∂Λ/∂φ(equation 3)
 end					

i = 1;
 for j = 1:(numel(th2(i, :)))
  N_c2 = numel(th2(i,j))/th2(i,j)+N_c2;		% "Stress Lv. 2" in the first term of ∂Λ/∂φ(equation 3)
 end	

i = 1;
 for j = 1:(numel(th3(i, :)))
  N_c3 = numel(th3(i,j))/th3(i,j)+N_c3;		% "Stress Lv. 3" in the first term of ∂Λ/∂φ(equation 3)
 end	

N_c = N_c1 + N_c2 + N_c3;				  % ∂Λ/∂φ first term

three = -B*(N_c  - (c1 + c2 + c3));		  	
% vpa(three,3);

%------------- equation 3 -------------%
d1 = 0;
d2 = 0;
d3 = 0;

i = 1;					% "Stress Lv. 1" in the second term of ∂Λ/∂b(equation 4)
 for j = 1:(numel(th1(i, :)))
  N = numel(th1(i,j));			
  d1 = N/th1(i+1,j)*(th1(i+2,j)/A*exp(-(phi/th1(i,j)+b/th1(i+1,j))))^B+d1;
 end

i = 1;					% "Stress Lv. 2" in the second term of ∂Λ/∂b(equation 4)"
 for j = 1:(numel(th2(i, :)))
  N = numel(th2(i,j));			
  d2 = N/th2(i+1,j)*(th2(i+2,j)/A*exp(-(phi/th2(i,j)+b/th2(i+1,j))))^B+d2;
 end

i = 1;					% "Stress Lv. 3" in the second term of ∂Λ/∂b(equation 4)
 for j = 1:(numel(th3(i, :)))
  N = numel(th3(i,j));			
  d3 = N/th3(i+1,j)*(th3(i+2,j)/A*exp(-(phi/th3(i,j)+b/th3(i+1,j))))^B+d3;
 end


 
N_d1 = 0;				
N_d2 = 0;
N_d3 = 0;

i = 1;
 for j = 1:(numel(th1(i, :)))
  N_d1 = numel(th1(i,j))/th1(i+1,j)+N_d1;		% "Stress Lv. 1" in the first term of ∂Λ/∂b(equation 4)
 end					

i = 1;
 for j = 1:(numel(th2(i, :)))
  N_d2 = numel(th2(i,j))/th2(i+1,j)+N_d2;		% "Stress Lv. 2" in the first term of ∂Λ/∂b(equation 4)
 end	

i = 1;
 for j = 1:(numel(th3(i, :)))
  N_d3 = numel(th3(i,j))/th3(i+1,j)+N_d3;		% "Stress Lv. 3" in the first term of ∂Λ/∂b(equation 4)
 end	

N_d = N_d1 + N_d2 + N_d3;				  % ∂Λ/∂b first term

four = -B*(N_d - (d1 + d2 + d3));		  
% vpa(four,3);

one

one = 

two

two = 

●

three

three = 

four

four = 

%[B, A, phi, b] = solve(one,two,three,four,'Real',true)	 
F = one^2 + two^2 + three^2 + four^2;
G = matlabFunction(F, 'vars', {[A B b phi]})

G = function_handle with value:
    @(in1)((((exp(in1(:,3).*(-5.0./2.0)-in1(:,4)./3.98e+2).*1.08e+2)./in1(:,1)).^in1(:,2)-1.0).*(in1(:,3).*(5.0./2.0)+in1(:,4)./3.98e+2-log(1.08e+2./in1(:,1)))+(((exp(in1(:,3).*(-5.0./2.0)-in1(:,4)./3.98e+2).*1.23e+2)./in1(:,1)).^in1(:,2)-1.0).*(in1(:,3).*(5.0./2.0)+in1(:,4)./3.98e+2-log(1.23e+2./in1(:,1)))+(((exp(in1(:,3).*(-5.0./2.0)-in1(:,4)./3.98e+2).*1.66e+2)./in1(:,1)).^in1(:,2)-1.0).*(in1(:,3).*(5.0./2.0)+in1(:,4)./3.98e+2-log(1.66e+2./in1(:,1)))+(((exp(in1(:,3).*(-5.0./4.0)-in1(:,4)./3.78e+2).*1.9e+2)./in1(:,1)).^in1(:,2)-1.0).*(in1(:,3).*(5.0./4.0)+in1(:,4)./3.78e+2-log(1.9e+2./in1(:,1)))+(((exp(in1(:,3).*(-5.0./4.0)-in1(:,4)./3.78e+2).*2.08e+2)./in1(:,1)).^in1(:,2)-1.0).*(in1(:,3).*(5.0./4.0)+in1(:,4)./3.78e+2-log(2.08e+2./in1(:,1)))+(((exp(in1(:,3).*(-5.0./2.0)-in1(:,4)./3.98e+2).*2.0e+2)./in1(:,1)).^in1(:,2)-1.0).*(in1(:,3).*(5.0./2.0)+in1(:,4)./3.98e+2-log(2.0e+2./in1(:,1)))+(((exp(in1(:,3).*(-5.0./4.0)-in1(:,4)./3.78e+2).*2.3e+2)./in1(:,1)).^in1(:,2)-1.0).*(in1(:,3).*(5.0./4.0)+in1(:,4)./3.78e+2-log(2.3e+2./in1(:,1)))+(((exp(in1(:,3).*(-5.0./4.0)-in1(:,4)./3.78e+2).*2.98e+2)./in1(:,1)).^in1(:,2)-1.0).*(in1(:,3).*(5.0./4.0)+in1(:,4)./3.78e+2-log(2.98e+2./in1(:,1)))+(((exp(in1(:,3).*(-5.0./2.0)-in1(:,4)./3.78e+2).*3.1e+2)./in1(:,1)).^in1(:,2)-1.0).*(in1(:,3).*(5.0./2.0)+in1(:,4)./3.78e+2-log(3.1e+2./in1(:,1)))+(((exp(in1(:,3).*(-5.0./2.0)-in1(:,4)./3.78e+2).*3.16e+2)./in1(:,1)).^in1(:,2)-1.0).*(in1(:,3).*(5.0./2.0)+in1(:,4)./3.78e+2-log(3.16e+2./in1(:,1)))+(((exp(in1(:,3).*(-5.0./2.0)-in1(:,4)./3.78e+2).*3.29e+2)./in1(:,1)).^in1(:,2)-1.0).*(in1(:,3).*(5.0./2.0)+in1(:,4)./3.78e+2-log(3.29e+2./in1(:,1)))+(((exp(in1(:,3).*(-5.0./2.0)-in1(:,4)./3.78e+2).*4.11e+2)./in1(:,1)).^in1(:,2)-1.0).*(in1(:,3).*(5.0./2.0)+in1(:,4)./3.78e+2-log(4.11e+2./in1(:,1)))+1.2e+1./in1(:,2)).^2+in1(:,2).^2.*(((exp(in1(:,3).*(-5.0./2.0)-in1(:,4)./3.98e+2).*1.08e+2)./in1(:,1)).^in1(:,2).*(5.0./2.0)+((exp(in1(:,3).*(-5.0./2.0)-in1(:,4)./3.98e+2).*1.23e+2)./in1(:,1)).^in1(:,2).*(5.0./2.0)+((exp(in1(:,3).*(-5.0./2.0)-in1(:,4)./3.98e+2).*1.66e+2)./in1(:,1)).^in1(:,2).*(5.0./2.0)+((exp(in1(:,3).*(-5.0./4.0)-in1(:,4)./3.78e+2).*1.9e+2)./in1(:,1)).^in1(:,2).*(5.0./4.0)+((exp(in1(:,3).*(-5.0./4.0)-in1(:,4)./3.78e+2).*2.08e+2)./in1(:,1)).^in1(:,2).*(5.0./4.0)+((exp(in1(:,3).*(-5.0./2.0)-in1(:,4)./3.98e+2).*2.0e+2)./in1(:,1)).^in1(:,2).*(5.0./2.0)+((exp(in1(:,3).*(-5.0./4.0)-in1(:,4)./3.78e+2).*2.3e+2)./in1(:,1)).^in1(:,2).*(5.0./4.0)+((exp(in1(:,3).*(-5.0./4.0)-in1(:,4)./3.78e+2).*2.98e+2)./in1(:,1)).^in1(:,2).*(5.0./4.0)+((exp(in1(:,3).*(-5.0./2.0)-in1(:,4)./3.78e+2).*3.1e+2)./in1(:,1)).^in1(:,2).*(5.0./2.0)+((exp(in1(:,3).*(-5.0./2.0)-in1(:,4)./3.78e+2).*3.16e+2)./in1(:,1)).^in1(:,2).*(5.0./2.0)+((exp(in1(:,3).*(-5.0./2.0)-in1(:,4)./3.78e+2).*3.29e+2)./in1(:,1)).^in1(:,2).*(5.0./2.0)+((exp(in1(:,3).*(-5.0./2.0)-in1(:,4)./3.78e+2).*4.11e+2)./in1(:,1)).^in1(:,2).*(5.0./2.0)-2.5e+1).^2+in1(:,2).^2.*(((exp(in1(:,3).*(-5.0./2.0)-in1(:,4)./3.98e+2).*1.08e+2)./in1(:,1)).^in1(:,2)./3.98e+2+((exp(in1(:,3).*(-5.0./2.0)-in1(:,4)./3.98e+2).*1.23e+2)./in1(:,1)).^in1(:,2)./3.98e+2+((exp(in1(:,3).*(-5.0./2.0)-in1(:,4)./3.98e+2).*1.66e+2)./in1(:,1)).^in1(:,2)./3.98e+2+((exp(in1(:,3).*(-5.0./4.0)-in1(:,4)./3.78e+2).*1.9e+2)./in1(:,1)).^in1(:,2)./3.78e+2+((exp(in1(:,3).*(-5.0./4.0)-in1(:,4)./3.78e+2).*2.08e+2)./in1(:,1)).^in1(:,2)./3.78e+2+((exp(in1(:,3).*(-5.0./2.0)-in1(:,4)./3.98e+2).*2.0e+2)./in1(:,1)).^in1(:,2)./3.98e+2+((exp(in1(:,3).*(-5.0./4.0)-in1(:,4)./3.78e+2).*2.3e+2)./in1(:,1)).^in1(:,2)./3.78e+2+((exp(in1(:,3).*(-5.0./4.0)-in1(:,4)./3.78e+2).*2.98e+2)./in1(:,1)).^in1(:,2)./3.78e+2+((exp(in1(:,3).*(-5.0./2.0)-in1(:,4)./3.78e+2).*3.1e+2)./in1(:,1)).^in1(:,2)./3.78e+2+((exp(in1(:,3).*(-5.0./2.0)-in1(:,4)./3.78e+2).*3.16e+2)./in1(:,1)).^in1(:,2)./3.78e+2+((exp(in1(:,3).*(-5.0./2.0)-in1(:,4)./3.78e+2).*3.29e+2)./in1(:,1)).^in1(:,2)./3.78e+2+((exp(in1(:,3).*(-5.0./2.0)-in1(:,4)./3.78e+2).*4.11e+2)./in1(:,1)).^in1(:,2)./3.78e+2-3.121427242030257e-2).^2+1.0./in1(:,1).^2.*in1(:,2).^2.*(((exp(in1(:,3).*(-5.0./2.0)-in1(:,4)./3.98e+2).*1.08e+2)./in1(:,1)).^in1(:,2)+((exp(in1(:,3).*(-5.0./2.0)-in1(:,4)./3.98e+2).*1.23e+2)./in1(:,1)).^in1(:,2)+((exp(in1(:,3).*(-5.0./2.0)-in1(:,4)./3.98e+2).*1.66e+2)./in1(:,1)).^in1(:,2)+((exp(in1(:,3).*(-5.0./4.0)-in1(:,4)./3.78e+2).*1.9e+2)./in1(:,1)).^in1(:,2)+((exp(in1(:,3).*(-5.0./4.0)-in1(:,4)./3.78e+2).*2.08e+2)./in1(:,1)).^in1(:,2)+((exp(in1(:,3).*(-5.0./2.0)-in1(:,4)./3.98e+2).*2.0e+2)./in1(:,1)).^in1(:,2)+((exp(in1(:,3).*(-5.0./4.0)-in1(:,4)./3.78e+2).*2.3e+2)./in1(:,1)).^in1(:,2)+((exp(in1(:,3).*(-5.0./4.0)-in1(:,4)./3.78e+2).*2.98e+2)./in1(:,1)).^in1(:,2)+((exp(in1(:,3).*(-5.0./2.0)-in1(:,4)./3.78e+2).*3.1e+2)./in1(:,1)).^in1(:,2)+((exp(in1(:,3).*(-5.0./2.0)-in1(:,4)./3.78e+2).*3.16e+2)./in1(:,1)).^in1(:,2)+((exp(in1(:,3).*(-5.0./2.0)-in1(:,4)./3.78e+2).*3.29e+2)./in1(:,1)).^in1(:,2)+((exp(in1(:,3).*(-5.0./2.0)-in1(:,4)./3.78e+2).*4.11e+2)./in1(:,1)).^in1(:,2)+1.2e+1).^2

format long g
[result, fval] = fmincon(G, [.1 .2 .3 .4])

Local minimum possible. Constraints satisfied.

fmincon stopped because the size of the current step is less than
the value of the step size tolerance and constraints are 
satisfied to within the value of the constraint tolerance.

result = 1×4

1.0e+00 * 223.521560318689 -0.128872549985852 18.6063254440608 -17463.55217938

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

fval = 
      0.000163646733037038

Those are the positions. The fval is not especially small, leading us to suspect that the minima might be greater than zero, indicating that potentially there is not solution to the problem.

G12 = @(x1,x2) G([x1, x2, result(3:4)]);
fsurf(G12, [220 240 -0.1289 -0.1288])

Warning: Function behaves unexpectedly on array inputs. To improve performance, properly vectorize your function to return an output with the same size and shape as the input arguments.

figure
G34 = @(x3,x4) G([result(1:2), x3, x4]);
[X3, X4] = ndgrid(linspace(18.58, 18.62,20), linspace(-17500,-17400,21));
surf(X3, X4, arrayfun(G34, X3, X4))

These graphs show that the location does seem to be a minima, and it doesn't look as-if the minima is ever 0 (0 would indicate a perfect fit for the four equations.)

[J]

Thank you both for taking the time to answer. It is judged to be a convergence problem. So is the partial differential of the equation wrong? Or is it the coding wrong? I was wondering if you could suggest a solution for each.

Answer 1

[Torsten]

If Λ is some kind of maximum likelihood function you try to maximize, I suggest you use "mle".

We had this problem already - I think there is nothing new that could be said if you can't give an underlying probability density function that leads to Λ:

https://uk.mathworks.com/matlabcentral/answers/2015586-error-on-symbolic-calculation-empty-sym-0-by-1?s_tid=prof_contriblnk

But your estimated parameters don't seem to solve the four equations:

B0 = 5.874395;
A0 = 6e-5;
phi0 = 5.630329851e3;
b0 = 2.80599e-1;

p0 = [B0;A0;phi0;b0];
norm(fun(p0))

ans = 2.3158e+06

format long
p = fsolve(@fun,p0,optimset('MaxFunEvals',1000000,'MaxIter',1000000))

No solution found.

fsolve stopped because the relative size of the current step is less than the
value of the step size tolerance squared, but the vector of function values
is not near zero as measured by the value of the function tolerance.

p = 4×1

1.0e+05 * 0.000000039232845 0.000016860827754 4.993994189282574 -0.004703145086915

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

norm(fun(p))

ans = 
   0.060744909668088

function res = fun(p)
B = p(1);
A = p(2);
phi = p(3);
b = p(4);

th1 = [378 378 378 378; 0.4 0.4 0.4 0.4; 310 316 329 411]; % [V_i ; U_i ; T_i] 
th2 = [378 378 378 378; 0.8 0.8 0.8 0.8; 190 208 230 298];
th3 = [398 398 398 398; 0.4 0.4 0.4 0.4; 108 123 166 200];
%syms B;
%syms A;
%syms phi;
%syms b;
%------------- equation 1 -------------%
a1_th1 = 0;					
a1_th2 = 0;
a1_th3 = 0;
i = 1;					% "Stress Lv. 1" in the second term of ∂Λ/∂β(equation 1)
 for j = 1:(numel(th1(i, :)))
  N = numel(th1(i,j));			
                                                	
  a1_th1 = N*(1-(th1(i+2,j)/A*exp(-(phi/th1(i,j)+b/th1(i+1,j))))^B)*(log(th1(i+2,j)/A)-(phi/th1(i,j)+b/th1(i+1,j))) + a1_th1 ;
	
 end
i = 1;					% "Stress Lv. 2" in the second term of ∂Λ/∂β(equation 1)
 for j = 1:(numel(th2(i, :)))
  N = numel(th2(i,j));			
                                                	
  a1_th2 = N*(1-(th2(i+2,j)/A*exp(-(phi/th2(i,j)+b/th2(i+1,j))))^B)*(log(th2(i+2,j)/A)-(phi/th2(i,j)+b/th2(i+1,j))) + a1_th2 ;
	
 end
i = 1;					% "Stress Lv. 3" in the second term of ∂Λ/∂β(equation 1)
 for j = 1:(numel(th3(i, :)))
  N = numel(th3(i,j));			
                                                	
  a1_th3 = N*(1-(th3(i+2,j)/A*exp(-(phi/th3(i,j)+b/th3(i+1,j))))^B)*(log(th3(i+2,j)/A)-(phi/th3(i,j)+b/th3(i+1,j))) + a1_th3 ;
	
 end
a1 = a1_th1 + a1_th2 + a1_th3;
N = numel(th1(1,:)) + numel(th2(1,:)) + numel(th3(1,:));		% ∂Λ/∂β first term
one = N/B + a1;			
% vpa(one,3);
%------------- equation 2 -------------%
b1 = 0;
b2 = 0;
b3 = 0;
i = 1;					% "Stress Lv. 1" in the second term of ∂Λ/∂A(equation 2)
 for j = 1:(numel(th1(i, :)))
  N = numel(th1(i,j));
  b1 = N*(th1(i+2,j)/A*exp(-(phi/th1(i,j)+b/th1(i+1,j))))^B+b1;	
 end
i = 1;					% "Stress Lv. 2" in the second term of ∂Λ/∂A(equation 2)
 for j = 1:(numel(th2(i, :)))
  N = numel(th2(i,j));
  b2 = N*(th2(i+2,j)/A*exp(-(phi/th2(i,j)+b/th2(i+1,j))))^B+b2;	
 end
i = 1;					% "Stress Lv. 3" in the second term of ∂Λ/∂A(equation 2)
 for j = 1:(numel(th3(i, :)))
  N = numel(th3(i,j));
  b3 = N*(th3(i+2,j)/A*exp(-(phi/th3(i,j)+b/th3(i+1,j))))^B+b3;	
 end
N = numel(th1(1,:)) + numel(th2(1,:)) + numel(th3(1,:));		% ∂Λ/∂A first term
two = -B/A*(N + (b1 + b2 + b3));		  	
% vpa(two,3);
%------------- equation 3 -------------%
c1 = 0;
c2 = 0;
c3 = 0;
i = 1;					% "Stress Lv. 1" in the second term of ∂Λ/∂φ(equation 3)
 for j = 1:(numel(th1(i, :)))
  N = numel(th1(i,j));			
  c1 = N/th1(i,j)*(th1(i+2,j)/A*exp(-(phi/th1(i,j)+b/th1(i+1,j))))^B+c1;
 end
i = 1;					% "Stress Lv. 2" in the second term of ∂Λ/∂φ(equation 3)
 for j = 1:(numel(th2(i, :)))
  N = numel(th2(i,j));			
  c2 = N/th2(i,j)*(th2(i+2,j)/A*exp(-(phi/th2(i,j)+b/th2(i+1,j))))^B+c2;
 end
i = 1;					% "Stress Lv. 3" in the second term of ∂Λ/∂φ(equation 3)
 for j = 1:(numel(th3(i, :)))
  N = numel(th3(i,j));			
  c3 = N/th3(i,j)*(th3(i+2,j)/A*exp(-(phi/th3(i,j)+b/th3(i+1,j))))^B+c3;
 end
N_c1 = 0;				
N_c2 = 0;
N_c3 = 0;
i = 1;
 for j = 1:(numel(th1(i, :)))
  N_c1 = numel(th1(i,j))/th1(i,j)+N_c1;		% "Stress Lv. 1" in the first term of ∂Λ/∂φ(equation 3)
 end					
i = 1;
 for j = 1:(numel(th2(i, :)))
  N_c2 = numel(th2(i,j))/th2(i,j)+N_c2;		% "Stress Lv. 2" in the first term of ∂Λ/∂φ(equation 3)
 end	
i = 1;
 for j = 1:(numel(th3(i, :)))
  N_c3 = numel(th3(i,j))/th3(i,j)+N_c3;		% "Stress Lv. 3" in the first term of ∂Λ/∂φ(equation 3)
 end	
N_c = N_c1 + N_c2 + N_c3;				  % ∂Λ/∂φ first term
three = -B*(N_c  - (c1 + c2 + c3));		  	
% vpa(three,3);
%------------- equation 3 -------------%
d1 = 0;
d2 = 0;
d3 = 0;
i = 1;					% "Stress Lv. 1" in the second term of ∂Λ/∂b(equation 4)
 for j = 1:(numel(th1(i, :)))
  N = numel(th1(i,j));			
  d1 = N/th1(i+1,j)*(th1(i+2,j)/A*exp(-(phi/th1(i,j)+b/th1(i+1,j))))^B+d1;
 end
i = 1;					% "Stress Lv. 2" in the second term of ∂Λ/∂b(equation 4)"
 for j = 1:(numel(th2(i, :)))
  N = numel(th2(i,j));			
  d2 = N/th2(i+1,j)*(th2(i+2,j)/A*exp(-(phi/th2(i,j)+b/th2(i+1,j))))^B+d2;
 end
i = 1;					% "Stress Lv. 3" in the second term of ∂Λ/∂b(equation 4)
 for j = 1:(numel(th3(i, :)))
  N = numel(th3(i,j));			
  d3 = N/th3(i+1,j)*(th3(i+2,j)/A*exp(-(phi/th3(i,j)+b/th3(i+1,j))))^B+d3;
 end
N_d1 = 0;				
N_d2 = 0;
N_d3 = 0;
i = 1;
 for j = 1:(numel(th1(i, :)))
  N_d1 = numel(th1(i,j))/th1(i+1,j)+N_d1;		% "Stress Lv. 1" in the first term of ∂Λ/∂b(equation 4)
 end					
i = 1;
 for j = 1:(numel(th2(i, :)))
  N_d2 = numel(th2(i,j))/th2(i+1,j)+N_d2;		% "Stress Lv. 2" in the first term of ∂Λ/∂b(equation 4)
 end	
i = 1;
 for j = 1:(numel(th3(i, :)))
  N_d3 = numel(th3(i,j))/th3(i+1,j)+N_d3;		% "Stress Lv. 3" in the first term of ∂Λ/∂b(equation 4)
 end	
N_d = N_d1 + N_d2 + N_d3;				  % ∂Λ/∂b first term
four = -B*(N_d - (d1 + d2 + d3));		  
% vpa(four,3);
res = [one;two;three;four];
end

[J]

The link you gave me is also my old question. You provided a smart solution. Thank you.

It's a nonlinear problem, but I wonder if I need to present approximated solution.

As long as there is no nonlinear function in the denominator of the partial differential equation, I think it can be solved with a linear system of equations, I would like to ask if you can think a little more.

[Torsten]

As long as there is no nonlinear function in the denominator of the partial differential equation, I think it can be solved with a linear system of equations.

I don't know what you mean here.

Maybe you should try what you get if you define only Λ as a function of the parameters and use "fmincon" to maximize it.

[J]

This means that there are no exponential functions, power functions, or logarithmic functions in the denominator.

Considering that you tipped me before, you told me to use mle, but the pdf of this likelihood function is as follows. I was wondering if you could give me any advice on how to apply mle here.

[Torsten]

There are two different kinds of parameter types that your T-H Weibull depends on, namely external parameters (V,U), distribution parameters (beta,A,phi,b) (to be fitted) and the observed quantity t. I don't see a way to use "mle" for distributions that are additionally parametrized by external (given) parameters.

So I think the way you proceeded is the only possible: explicitly write down the log-likelihood function Λ for your data. But after that, I would first try to apply "fmincon" to maximize Λ with beta,A,phi and b being the parameters to be fitted. If this does not work, do as you already tried: differentiate the log-likelihood function with respect to the parameters, set the expressions to zero and try to solve the four nonlinear equations for the parameters.

Obviously, you made a mistake in your partial derivatives - the following code gives your supplied solution β = 5.874395, A = 0.000060, b = 0.280599, φ = 5630.329851. Note that extrema for f are also extrema for log(f). You may try maximizing log(f) instead of f as I did below, but it took too long to see whether this computation was successful, too.

th1 = [378 378 378 378; 0.4 0.4 0.4 0.4; 310 316 329 411]; % [V_i ; U_i ; T_i] 
th2 = [378 378 378 378; 0.8 0.8 0.8 0.8; 190 208 230 298];
th3 = [398 398 398 398; 0.4 0.4 0.4 0.4; 108 123 166 200];
th = [th1,th2,th3];
V = th(1,:);
U = th(2,:);
T = th(3,:);
syms beta A b phi
f = prod(beta/A*exp(-(phi./V+b./U)).*(T/A.*exp(-(phi./V+b./U))).^(beta-1).*exp(-(T/A.*exp(-(phi./V+b./U))).^beta));
dfdbeta = diff(f,beta);
dfdA = diff(f,A);
dfdb = diff(f,b);
dfdphi = diff(f,phi);
format long
[beta_sol,A_sol,b_sol,phi_sol] = vpasolve([dfdbeta,dfdA,dfdb,dfdphi]==[0,0,0,0],[beta,A,b,phi],[5.87, 0.00006, 0.28, 5630.3])

beta_sol = 

5.8744444936953754728727701062122

A_sol = 

0.000059702012870644238396009145567829

b_sol = 

0.28059831624318766178811912460829

phi_sol = 

5630.3264041582504081439236572513

double(subs([dfdbeta,dfdA,dfdb,dfdphi],[beta,A,b,phi],[beta_sol,A_sol,b_sol,phi_sol]))

ans = 1×4

1.0e-60 * 0.000001483682460 -0.933452291679171 0.000151929083932 -0.000000148368246

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

●