find the root of a function in an interval

Hi fellows,

I would like to solve an expontential function f in an interval, say, [0,e]. I try to use solve but it is very slow. Then I try to use fzero .It is faster but still not efficient enough. As I want the solution to be in the interval, I think it should be more efficient if I can specifiy the interval. Luckily fzeros allows me to do so. But meanwhile it requires the value of function f changes sign in the boundadry of the interval. I can not guarantee f(0) has a different sign as f(e) so fzero is not working...Could you give me some advice about this? Here is an example of the code

if true

    % 
r=0.01;
D=[0.1,-0.5,-0.4;-0.6,0.2,0.5];
P=[0.2,0.2,0.6];
d11=D(1,1);
d12=D(1,2);
d13=D(1,3);
p1=P(1);
p2=P(2);
p3=P(3);
e=100;
f=@(x) p1*d11*exp(-r*x*d11)+p2*d12*exp(-r*x*d12)+p3*d13*exp(-r*x*d13);
x=fzero(f,[0,100])
  end

and I get the error messages is

if true
  % Error using fzero (line 274)
The function values at the interval endpoints must differ in sign.

end

[Andrei Bobrov]

Please show your function

[xueqi]

Yes I just did it.

[dpb]

Well, the only way it can be <0 is for one or more of the P be negative and of absolute value >(sum of other two). But since you've given them fixed values >0 there is no zero possible. exp(-z)-->0 as z-->inf but it's not actually zero until you reach underflow on the smallest of the three terms.

That you can actually calculate from the limits of FP precision.

The function is 0.2*{1.1->0} + 0.2*{0.5->0} + 0.6*{0.6->0} as a range of x

ADDENDUM:

Given your above values,

>> D=[1.1,0.5,0.6;0.4,1.2,1.5];
>> P=[0.2,0.2,0.6];
>> r=D(1,:).*P;
>> m=eps(0)./r
m =
9.9e-323 *
  0.2500    0.5000    0.1500
>> log(m)
ans =
-742.8306 -742.1375 -743.3415
>>

>> The value of the argument of the exponential for each of the three terms can be no larger (in absolute value) than the three values above and avoid underflow. So, you can compute the point of the total expression going to look like zero as those values where rx*D1j is the corresponding value.

ADDENDUM 2:

I'm guessing the above isn't exactly what you intended; it simply reflects the algebra of what you posted. Your properly-formulated problem is likely quite different... :)

[xueqi]

Hi, Thanks for your reply. I actually did make a mistake specifying D. Now it is possible for f to be negative. But what I really want to is find a general way to do this. As I have a set of D that need to be calculated, I am unable to analyze them individually.

[dpb]

A set of "D" that need to be calculated???? You've given D as a set of constants.

You need to sit down with pencil and paper and look at what you get for limiting conditions from your set of constants -- or plot it and visualize it. The system is not at all like you think methinks...

Try

>> P=[0.2,0.2,0.6];
>> D=[0.1,-0.5,-0.4];
>> r=0.01;
>> rD=r*D;
>> f=@(x) P.*D.*exp(-rD.*x);
>> x=[-100:100]';
>> y=cell2mat(arrayfun(f,x,'uniformoutput',false));
>> plot(x,[y sum(y,2)])
>>