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!
采纳的回答
Walter Roberson
2017-2-14
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
7 个评论
Egle
2017-2-14
Egle
2017-2-18
Walter Roberson
2017-2-18
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
2017-3-16
Walter Roberson
2017-3-16
Your function declaration has no slot to receive D or R
Egle
2017-3-17
Walter Roberson
2017-3-17
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 );
更多回答（2 个）
Torsten
2017-3-17
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.
3 个评论
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.
Torsten
2017-3-22
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.
类别
在 帮助中心 和 File Exchange 中查找有关 Loops and Conditional Statements 的更多信息
产品
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
