[Scilab-users] root calculation

[Calixte Denizet]

Hi Paul,

If the exponent are always negative and the coefficients are always 
positive, the function is strictly decreasing (so the derivative is 
non-null), you can use a Newton method.

function y=newton(fun, dfun, x0, eps)
     y=fun(x0)
     while abs(y) > eps
         x0 = x0 - y / dfun(x0);
         y = fun(x0)
     end
     y = x0;
endfunction

deff('y=foo(x)','y=0.403*x.^(-0.121)+60.5*x.^(-0.73)-0.1839')
deff('y=dfoo(x)','y=0.403*-0.121*x.^(-1.121)+60.5*-0.73*x.^(-1.73)')

To find a "good" starting point:
If x is a solution then 0.403*x^(-0.121)<=0.1839 and 
60.5*x^(-0.73)<=0.1839, so x >=max((0.1839/0.403)^(-1/0.121), 
(0.1839/60.5)^(-1/0.73))

newton(foo,dfoo,2806,1e-10)

A faster way:
function y=newton2(coeffs, expo, cste, x0, eps)
     y = coeffs * (x0 .^ expo)' - cste;
     dcoeffs = coeffs .* expo;
     dexpo = expo - 1;
     while abs(y) > eps
         x0 = x0 - y / (dcoeffs * (x0 .^ dexpo)');
         y = coeffs * (x0 .^ expo)' - cste;
     end
     y = x0;
endfunction

newton2([0.403 60.5], [-0.121 -0.73],0.1839,2806,1e-10)

Calixte

On 05/11/2012 13:13, Paul Carrico wrote:
>
> Dear all,
>
> This a stupid question, but how can I solve directly in, Scilab  an 
> equation such as :
>
> 0.403*X^(-0.121) + 60.5*X^(-0.73) -- 0.1839 = 0
>
> ?
>
> Is-it necessary to code a function ? from memory : dichotomy method, 
> secant method, Brent one etc. ...
>
> Regards
>
> Paul
>
>
>
> _______________________________________________
> users mailing list
> users at lists.scilab.org
> http://lists.scilab.org/mailman/listinfo/users


-- 
Calixte Denizet
Software Development Engineer
-----------------------------------------------------------
Scilab Enterprises
143bis rue Yves Le Coz - 78000 Versailles, France
http://www.scilab-enterprises.com

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.scilab.org/pipermail/users/attachments/20121105/68ee3010/attachment.htm>