solving equation from vpasolve or anyway of solving

Hello, I have an oscillatory decaying equation f(h), when it intersects with zero, it should give me multiple values, but vpasolve or fsolve is only giving one value close to the initial value I provided, but I prefer to get all the solution in the range from [25, 50]. The equation is below: syms h P0=4.181323907591596e+03; d1=13.5; d=9.734259414551886; d2=4.706483519517566; F=(2*pi*P0*sin(2*pi*(h-d1)/d1)*exp(d^3/(d1^2*d2)-(h-d1)/d2)/(d1*10e-9)+... P0*cos(2*pi*(h-d1)/d1)*exp(d^3/(d1^2*d2)-(h-d1)/d2)/(d2*10e-9))==0; S=vpasolve(F,h,[25,50])

If you can tell me where I went wrong, that would be great, thanks. PS: You can try to change h from [10,50] to see the intersected values, but I am not able to obtain each one of them numerically.

Answer 1

[Stephan]

Hi,

you can specify the interval for h which is used for the solution:

syms h 
P0=4.181323907591596e+03; 
d1=13.5;
d=9.734259414551886; 
d2=4.706483519517566;
F=(2*pi*P0*sin(2*pi*(h-d1)/d1)*exp(d^3/(d1^2*d2)-(h-d1)/d2)/(d1*10e-9)+...
    P0*cos(2*pi*(h-d1)/d1)*exp(d^3/(d1^2*d2)-(h-d1)/d2)/(d2*10e-9))==0;
sol = (solve([F, h >= 25, h <= 50], h));
sol_num = double(sol)

This gives:

sol_num =

     32.8298
     26.0798
     46.3298
     39.5798

To check this results you can insert the results in a function handle of your function and expect that it will result in F=0 (note that this is calculated with the symbolic results sol - not sol_num:

fun = @(P0, h, d, d1, d2)(2*pi*P0.*sin(2*pi.*(h-d1)./d1).*exp(d^3/(d1^2*d2)-(h-d1)/d2)/(d1*10e-9)+...
    P0.*cos(2*pi*(h-d1)/d1).*exp(d^3/(d1^2*d2)-(h-d1)/d2)./(d2*10e-9))
result_numeric = double(fun(P0,sol,d,d1,d2))

The results are:

result_numeric =

     1.0e-61 *

      0.0668
      0.3160
      0.0198
      0.0061

which is pretty near by 0 and meets the conditions above. If you do the same calculation with the numeric results:

fun = @(P0, h, d, d1, d2)(2*pi*P0.*sin(2*pi.*(h-d1)./d1).*exp(d^3/(d1^2*d2)-(h-d1)/d2)/(d1*10e-9)+...
    P0.*cos(2*pi*(h-d1)/d1).*exp(d^3/(d1^2*d2)-(h-d1)/d2)./(d2*10e-9))
result_numeric = double(fun(P0,sol_num,d,d1,d2))

the result is:

result_numeric =

     1.0e-04 *

     -0.1478
     -0.4768
     -0.0069
      0.0334

due to rounding errors - but still a good result.

Best regards

Stephan

[Subinuer]

Thanks Stephan!