Using fzero function to solve a nonlinear equation with one input

Hi,

I am trying to solve this non-linear equation for alpha using a range of values of X (0-100) [so that I can plot alpha vs X].

function y=equations(alpha)

[D, A, R] = geom();
[n, m] = setglbl();
[X] = setX();

              Ag            = (R.^2)*(pi-alpha+0.5*sin(2*alpha))             ; 
              Sgw           = 2*R*(pi-alpha)                                 ; 
              Si            = 2*R*sin(alpha)                                 ; 
              Dg            = 4*Ag/(Sgw+Si)                                  ; 
              Al            = (R.^2)*(alpha-0.5*sin(2*alpha))                ; 
              Slw           = 2*R*alpha                                      ; 
              Dl            = 4*Al/Slw                                       ; 
              ag            = Ag/A                                           ; 

y=((X.^2).*(1-ag).^(n-2)*(D/Dl).^(n+1))-(ag.^(m-2)*(D/Dg).^(m+1)*(1+(Dg*Si/Al)));

end

WHERE

function [X] = setX
global X ;
X=linspace(0,100,101) ;
end

AND

function [Cg, Cl, n, m] = setglbl
global Cg Cl n m
n=1;
m=1;
end

AND

function [D, A, R] = geom()
D = 0.3; %Pipe Diameter (m)
R = D/2; %Pipe Radius (m)
A = pi*(R^2); %Cross-sectional Area (m^2)
end

I KEEP GETTING ERRORS:

Error using fzero (line 301)
FZERO cannot continue because user-supplied function_handle ==> equations failed with the error below.

   *Matrix dimensions must agree.*
  >> equations

*Not enough input arguments.*

Error in equations (line 7)
            Ag            = (R.^2)*(pi-alpha+0.5*sin(2*alpha))             ; % [m^2] Cross-sectional area of gas phase in pipe

Does anyone have some advice on what the problem could be?

Should the X interval be defined differently? Or should X be used as a second input to the function?

Any help appreciated!

Answer 1

[Walter Roberson]

You need to change

function [n, m] = setglbl
global n m
n=1;
m=1;
end

However, you have

X=linspace(0,100,101) ;

so X is a vector of length 101, and you have

y=((X.^2).*(1-ag).^(n-2)*(D/Dl).^(n+1))-(ag.^(m-2)*(D/Dg).^(m+1)*(1+(Dg*Si/Al)));

so y is going to be a vector the same length as X. You cannot fzero an entire vector.

Your formula has some numeric problems that need to be worked around. Revised source is:

function alphaVsX
    [D, A, R] = geom();
    [n, m] = setglbl();
    X = setX();

      alpha = arrayfun( @(x) fzero( @(alpha) equations(alpha, x, D, A, R, m, n), [1E-4, pi-1E-4]), X );
      plot(X, alpha);
      xlabel('X')
      ylabel('alpha')
  end

function [D, A, R] = geom()
    D = 0.3; %Pipe Diameter (m)
    R = D/2; %Pipe Radius (m)
    A = pi*(R^2); %Cross-sectional Area (m^2)
end

function [n, m] = setglbl
    n=1;
    m=1;
end

function [X] = setX
    X=linspace(0,100,101) ;
end

function y = equations(alpha, X, D, A, R, m, n)
    Ag            = (R.^2)*(pi-alpha+0.5*sin(2*alpha))             ;
    Sgw           = 2*R*(pi-alpha)                                 ;
    Si            = 2*R*sin(alpha)                                 ;
    Dg            = 4*Ag/(Sgw+Si)                                  ;
    Al            = (R.^2)*(alpha-0.5*sin(2*alpha))                ;
    Slw           = 2*R*alpha                                      ;
    Dl            = 4*Al/Slw                                       ;
    ag            = Ag/A                                           ;

      if X == 0
          y = 0;
      else
          y=((X.^2).*(1-ag).^(n-2)*(D/Dl).^(n+1)) - (ag.^(m-2)*(D/Dg).^(m+1)*(1+(Dg*Si/Al)));
      end

end

[Egle]

This is extremely helpful, thank you Walter! I have spent weeks looking at different methods, but being new to Matlab, I don't think I would have come up with this. Thanks again!

[Egle]

Hi Walter,

Sorry to bother you again. Just wondering, how did you choose the interval [1E-4, pi-1E-4]?

I am trying to use this code for similar problems and it doesn't seem to work anymore.

Thanks! Egle

[Walter Roberson]

I studied the behavior of the function at values close to 0 and close to pi.

At close to 0, alpha-0.5*sin(2*alpha) (for Al) becomes 0, which makes Dl 0, which gives a problem because you divide by Dl.

At close to Pi, (pi-alpha+0.5*sin(2*alpha)) (for Ag) becomes 1, and then (1-ag) becomes 0 in y, and since n = 1, that makes (1-ag).^(n-2) becomes 0^(-1) or 1/0 which is a problem.

If you have other formula, you need to study their boundary conditions and figure out the places where there are infinities or divisions by 0 numerically rather than just algebraically.

For example, sin(alpha) near 0 becomes approximately alpha, so sin(2*alpha) becomes approximately 2*alpha, so 0.5*sin(2*alpha) becomes approximately 1/2*2*alpha = alpha, so alpha-0.5*sin(2*alpha) becomes approximately 0. As alpha approaches 0, there is going to be a point where this becomes true due to floating point round-off.

[Egle]

Hello again Walter,

I am trying to adjust previous code to solve a slightly different problem (with two inputs this time). But it does not seem to work.

Could you please give me some advice on this?

My aim is to solve a non-linear equation for alpha with given values of X and al (both ranges of values). Here is the code,

INPUTS:

function [n, m, Z] = setglblch
n=1   ;
m=1   ;
Z=14;
end

function [D, A, R] = geom()
D = 0.3; %Pipe Diameter (m)
R = D/2; %Pipe Radius (m)
A = pi*(R^2); %Cross-sectional Area (m^2)
end

function [X] = setX
X=linspace(0,100,101) ;
end

function [al] = setal
al=linspace(0,1,101) ;
end

SOLVING:

function y = equationsch(alpha, X, A, m, n, Z, al)

   beta1 = (1-al).^(-1) ;
   beta2 = al/((1/A)-al) ;
   beta3 = 1/alpha ;
   beta4 = beta1-beta2.*beta3 ;
   beta = beta4.^(-1) ;

      if X == 0
          y = 0;
      else
          y1 = al.^(1-0.5) ; 
          y2 = (1-al).^(0.5*m-1) ;
          y3 = X/Z ;
          y4 = beta.^(1+m) ;
          y5 = alpha.^(1+n) ; 
          y6 = (y4/y5).*y3 ;
          y = y1.*y2-y6 ;         
      end
  end

function [alpha] = getalphach
[A] = geom();
[n, m, Z] = setglblch();
X = setX();
al = setal();

    alpha = arrayfun( @(x) fzero( @(alpha) equationsch(alpha, X, D, A, R, m, n, Z, al), [1E-4, pi-1E-4]), X, al );
  end

I GET ERROR Too many input arguments.

Can this be fixed? Is it okay to use arrayfun in this case?

Thanks in advance!

Egle

[Walter Roberson]

Your function declaration has no slot to receive D or R

[Egle]

Hello again,

Thank you for your reply. Sorry, I didn't need D or R in the function, but it still does not work after deleting it. (Error: too many inputs)

Both my inputs X and al are vectors. Could this be causing the issue?

Thanks in advance! Egle

[Walter Roberson]

You have

      alpha = arrayfun( @(x) fzero( @(alpha) equationsch(alpha, X, A, m, n, Z, al), [1E-4, pi-1E-4]), X, al );

which has three arguments for arrayfun: the function handle, X, and al . arrayfun() takes a function handle and passes each of the remaining arguments to it, so because you are providing two arguments after the function handle, the function handle has to take two arguments. But it does not: the function handle takes only one argument, lower-case x. If each X entry is to be paired with each al entry, then you need to have your function handle

@(x) fzero( @(alpha) equationsch(alpha, X, A, m, n, Z, al), [1E-4, pi-1E-4])

take two arguments instead.

After that, notice that your function handle has a formal parameter lower-case x as an argument, but it does nothing with that argument. Instead, it passes upper-case X to the function, which is the entire vector of X values.

Perhaps you want

        alpha = arrayfun( @(one_X, one_al) fzero( @(alpha) equationsch(alpha, one_X, A, m, n, Z, one_al), [1E-4, pi-1E-4]), X, al );

Answer 2

[Torsten]

Does this work ?

function [alpha] = getalphach
[A] = geom();
[n, m, Z] = setglblch();
X = setX();
al = setal();

for i=1:numel(X)
  Xloc = X(i);
  for j = 1:numel(al)
    alloc = al(j);
    alpha(i,j) = fzero(@(x)equationsch(x,Xloc,A,m,n,Z,alloc),[1e-4, pi-1e-4]);
  end 
end

Best wishes

Torsten.

[Egle]

Hi Torsten,

Thank you for your help.

Now I get a different error:

Error using fzero (line 252) Function values at interval endpoints must be finite and real.

Error in getalphach (line 10) alpha(i,j) = fzero(@(x)equationsch(X,Xloc,m,n,Z,alloc),[1e-4, pi-1e-4]);

[Torsten]

Try to call fzero with a single starting value, e.g.

 alpha(i,j) = fzero(@(x)equationsch(x,Xloc,A,m,n,Z,alloc),pi-1e-4);

If fzero cannot find a solution, you should first plot your function to see whether a zero really exists.

By the way: Your call to fzero must be with a small x:

Not:

alpha(i,j) = fzero(@(x)equationsch(X,Xloc,m,n,Z,alloc),[1e-4, pi-1e-4]);

But:

alpha(i,j) = fzero(@(x)equationsch(x,Xloc,m,n,Z,alloc),[1e-4, pi-1e-4]);

Best wishes

Torsten.

[Egle]

Dear Torsten,

Thank you for your help! And apologies for disturbing you again.

I finally managed to get your suggested double loop working! The problem was that my al range was defined as: al=linspace(0,1,101) ; and the 0 was putting the code off, once changed to: al=linspace(0.01,1,101) ; it's working fine.

However, I now realised that I cannot use any combinations of X and al ranges. I have a set of experimental (coupled) values instead.

If al and X were defined as such,

function [X, al] = setmartvar

    datalockmart = [0.01 1 1.05 ; 0.02 1 1.07 ; 0.04 1 1.12 ; 0.07 0.96 1.19 ; 0.10 0.96 1.24 ; 0.2 0.91 1.4 ; 0.4 0.86 1.7 ; 0.7 0.81 2.16 ; 1.0 0.77 2.61 ; 2.0 0.69 4.12 ; 4.0 0.60 7 ; 7 0.52 11.2 ; 10 0.47 15 ; 20 0.34 27.3 ; 40 0.24 50 ; 70 0.16 82 ; 100 0.1 111] ;

    X = datalockmart(:,1) ;
    aglm = datalockmart(:,2) ;
    al=1-aglm ;
  end

Can this be incorporated in the loops instead?

function [alpha] = getalphach
[A] = geom();
[n, m, Z] = setglblch();
X = setX();
al = setal();

for i=1:numel(X)
Xloc = X(i);
  for j = 1:numel(al)
    alloc = al(j);
    alpha(i,j) = fzero(@(x)equationsch(x,Xloc,A,m,n,Z,alloc),[1e-4, pi-1e-4]);
  end 
end

Thanks in advance!

Egle

Answer 3

[Torsten]

Sure.

function [alpha] = getalphach
[A] = geom();
[n, m, Z] = setglblch();
[X,al] = setmartvar();
for i=1:numel(X)
  Xloc = X(i);
  alloc = al(i);
  alpha(i) = fzero(@(x)equationsch(x,Xloc,A,m,n,Z,alloc),[1e-4, pi-1e-4]);
end

Best wishes

Torsten.

[Egle]

This helps, thank you! The code is finally working as expected