Should I trust solutions obtained for a nonlinear non-transcendental equation in matlab using 'fzero' function?

y = @(x)((10^-9)*x*(cosh(x)*sin(x)-sinh(x)*cos(x)) - ((10^-9)-1*(x^4)*(1+cosh(x)*cos(x))))

Suppose I have a kind of equation similar to the one above:

'y' is a function in variable 'x', it is a nonlinear, non-transcendental eqaution.

I have tried 'fzero', and I hae also plotted the function. So that plot gives me a shot regarding approximate location of the point whereas the function crosses y =0 line (corresponding x point is the answer). This is how I found the initial guess (eyeballing) and found following values for x:

0.149524, 1.875236, 4.694101

But wenenver I apply these numbers into the main equaiton to find the value of function at the mentioned points using 'subs' function in matlab, it is NOT equal to zero!!!!! (however we know value of function at the solution points should be equal to zero!!!!). Of course, using funciton 'round' yields zero. In short, I am a bit hessitant if my roots are correct or not.

Does anybody havea recommendation to find the roots? (logically snice I plot the function, I can estimate the first root should be pretty close to 0.15, and applying this initial guess into the 'fzero' function, yields 0.149524.) I don't know should trust it or not. (another reasin which makes me hesitant is the fact that physically the system obtaining such kind of nonlinear non-transcendental eqaution should have the first root close to 1.87).

Does anybody has a recommendation to solve such an enigma?

Thanks a lot

[Walter Roberson]

[-10.99554073487622, -7.854757438239676, -4.694091132983881, -1.875104068844092, -0.004728708044671077, 0.004728708044671077, 1.875104068844092, 4.694091132983881, 7.854757438239676, 10.99554073487622]

[Walter Roberson]

Which release are you using? When I fzero(y,0.15) I get -1.87510406884409 not 0.149524

[Alireza Babaei]

Can you please run my code in the following:

clc
clear all
close all
syms x laL
L = 100e-3;
rm = input('rm = ')
rs = rm
x = linspace(0,10,100)
y1 = (rm*rs*(L^3)*laL*(cosh(laL)*sin(laL) - sinh(laL)*cos(laL)) - (rs*(L^3) - rm*(laL^4)*(1+cosh(laL)*cos(laL))))
y2 = (rm.*rs.*(L.^3).*x.*(cosh(x).*sin(x) - sinh(x).*cos(x)) - (rs.*(L^3) - rm.*(x.^4).*(1+cosh(x).*cos(x))))
xi = linspace(0,10,11)
for i = 1 : length(xi)
    aa = vpasolve(y1 == 0, xi(i))
end 
y3 = 1 + cos(x).*cosh(x)
subplot(2,1,1)
plot(x,y2,'LineWidth',2)
grid on
ylim([-.000001 .000001])
hold on l
plot(x,y3,'--r','LineWidth',2)
subplot(2,1,2)
plot(x,y2,'LineWidth',2)
grid on
ylim([-400 400])
hold on l
plot(x,y3,'--r','LineWidth',2)
hold off

Running this code, shows that the main function (presented in 'y1' or 'y2' has another root pretty much close to zero (almost 0.14)). I also plotted y function to see if there is really another root close to zero and it weems that it is. However practically it does not match with reality.

This is the thing making me confused.

On the other hand you say that your root even for initial guess around 0.15 is NOT 0.14

[Walter Roberson]

What value for rm?

Answer 1

[Matt J]

fzero is a numerical root finder. One doesn't expect y(x) to be exactly zero at the roots that it finds, but it should be pretty close. One way to evaluate the result is to plot the function in a tight interval around the supposed root and see how close to zero it is compared to its neighbors. The three roots you've indicated look pretty good to me.

fun = @(x)((10^-9)*x*(cosh(x)*sin(x)-sinh(x)*cos(x)) - ((10^-9)-1*(x^4)*(1+cosh(x)*cos(x))));
   
[x1,y1]=fzero(fun,[0,0.2])
[x2,y2]=fzero(fun,[1,2])
[x3,y3]=fzero(fun,[4.2,5])


close all 

figure;
fplot(fun, x1+[-1,1]*.0001);
hold on; plot(x1,y1,'ro','MarkerSize',6); hold off


figure;
fplot(fun, x2+[-1,1]*.0001);
hold on; plot(x2,y2,'ro','MarkerSize',6); hold off

figure;
fplot(fun, x3+[-1,1]*.0001);
hold on; plot(x3,y3,'ro','MarkerSize',6); hold off

[Alireza Babaei]

thanks for your answer,

I reached sth wierd, actual equation has coefficients that I introduced them as: 'input ('coefficient')',

I dont know why but whenever I assign the value of coefficient directly into the equaiton (just like the one I posted), it solved for 1.87, but whenenver I use 'inout' function, it gives 0.15

[Steven Lord]

For some functions, the difference between values that are displayed in the current display format and the values that are actually stored is negligible.

y = @(x) 2*x;
p = pi/2;
p2 = p + eps(p);
[p, y(p); p2, y(p2)]

In this case, p and p2 display the same (using the default format.) Their values are very close (but not the same, as you can see by computing p2-p) and the values computed by y are also very close.

y2 = @tan;
[p, y2(p); p2, y2(p2)]

A very small difference between p and p2 causes a huge difference in the value of y2. They even have different signs. Now this is a somewhat artificial example that I carefully concocted to have p and p2 end up on different sides of the tangent function's asymptote. But sometimes (especially when your function is steep) small differences get magnified.

If you plot some of the components of your function, when x gets large enough the cosh, sinh, and x.^4 functions start growing rapidly. So a very small change in x can result in a large change in y.

The fzero function is nice in that you can check if the answer is correct. Simply call the function you passed into fzero with the output of fzero as its input. If the answer is "close enough" to 0, the answer is a correct answer. If it's not, you may need to adjust the tolerances in your fzero call to define "close enough" how you want it to be defined.

[Matt J]

thanks for your answer

You're welcome. Please accept-click the answer if your question is resolved.

I dont know why but whenever I assign the value of coefficient directly into the equaiton (just like the one I posted), it solved for 1.87, but whenenver I use 'inout' function, it gives 0.15

You should plot both versions to see if the function is different in each case. Possibly your manual assignment specified the coefficient to more limited precision than the non-manual assignment.

[Walter Roberson]

0.149524 would be the solution if the 10^-9 in the equation (in both places) was 0.001

[Alireza Babaei]

In short, let me expose my experiece here:

I tried two different methods in finding the roots of the main function (y); using 'vpasolve' (and subsequently using plot to verify the results) and using 'fzero' (subsequently using fplot to verfy the obtianed roots).

Interesting point is: by means of vpasolve MATLAB shows a root pretty much close to the zero (0.14...). However using fzero, MATLAB does not show such a root and fplot does not depict/show any roots close to the zero neither.

This is means that plot verifies results obtained through vpasolve and/vs fplot verifies results shown by fzero (and these two (vpasolve and fzero) are incontrast to each other!)

If anyone finds another findings (othrthan whatever I mentioned) please let me know.

Now for me, it is a bit confusing which one to trust ...

[Matt J]

I suggest that you show us the plotting commands which reveal the presence of a root near 0.14. Below is my effort to plot the function this time using plot() with 100 million evenly spaced points in the interval shown. There is no indication of any root near 0.14.

fun = @(x)((10^-9)*x*(cosh(x)*sin(x)-sinh(x)*cos(x)) - ((10^-9)-1*(x^4)*(1+cosh(x)*cos(x))));
fun=str2func(vectorize(fun));

x=linspace(0.13,0.16,1e8);

plot(x,fun(x))

[Alireza Babaei]

Following is the plot I am talking about.

If you run the code I sent to you about 2 hours ago you'll get the following (just note that for rm I tried 0.001 (rm = 0.001))

BTW, first root (named as 'aa') is 0.14953642341834813206969708330693

[Matt J]

I don't see any code posted 2 hours ago. I suggest you re-post the code that generates the plots.

[Alireza Babaei]

sure here it is: run it and you can assign rm = 0.001 (for example)

clc
clear all
close all
syms x laL
L = 100e-3;
rm = input('rm = ')
rs = rm
x = linspace(0,10,100)
y1 = (rm*rs*(L^3)*laL*(cosh(laL)*sin(laL) - sinh(laL)*cos(laL)) - (rs*(L^3) - rm*(laL^4)*(1+cosh(laL)*cos(laL))))
y2 = (rm.*rs.*(L.^3).*x.*(cosh(x).*sin(x) - sinh(x).*cos(x)) - (rs.*(L^3) - rm.*(x.^4).*(1+cosh(x).*cos(x))))
xi = linspace(0,10,11)
for i = 1 : length(xi)
    aa = vpasolve(y1 == 0, xi(i))
end 
y3 = 1 + cos(x).*cosh(x)
subplot(2,1,1)
plot(x,y2,'LineWidth',2)
grid on
ylim([-.000001 .000001])
hold on l
plot(x,y3,'--r','LineWidth',2)
subplot(2,1,2)
plot(x,y2,'LineWidth',2)
grid on
ylim([-400 400])
hold on l
plot(x,y3,'--r','LineWidth',2)
hold off

[Matt J]

I reduced your code to the following.

L = 100e-3;
rm = 0.001; rs=rm;

y1 = @(laL) (rm.*rs.*(L.^3).*laL.*(cosh(laL).*sin(laL) - sinh(laL).*cos(laL)) ...
             - (rs.*(L.^3) - rm.*(laL.^4).*(1+cosh(laL).*cos(laL))));
         
         
fplot(y1,[0.13,0.16])         

[xr,yr]=fzero(y1,[0.13,0.16])

It does indeed have a root near 0.1495, according to both fplot and fzero. So I'm not sure what the issue is.

xr =

    0.1495


yr =

  -4.7069e-23

[Alireza Babaei]

Matt,

Yes, you're rigth.

I have no idea why this happens. Only thing comming to my mind refers to the coefficients of main function (which sounds ridiculus). maybe using input function yields the 0.14 but whenever applying the coefficient value directly into the main function, it does not show 0.14. (I am still not sure about it)

[Matt J]

It is almost certain that you are not giving the same values to the function in each case. Use a breakpoint to check the input values, or print them to the screen in long format, e.g.,

>> format long
>> pi

ans =

   3.141592653589793

[Walter Roberson]

syms laL
>> pretty(y1(laL))
   4
laL  (cos(laL) cosh(laL) + 1)   laL (cos(laL) sinh(laL) - cosh(laL) sin(laL))      1
----------------------------- - --------------------------------------------- - -------
             1000                                 1000000000                    1000000

>> pretty(y(laL))
   4                            laL (cos(laL) sinh(laL) - cosh(laL) sin(laL))        1
laL  (cos(laL) cosh(laL) + 1) - --------------------------------------------- - ----------
                                                  1000000000                    1000000000

They are not the same equation. The y one has first term multiplied by 1000 and the third term divided by 1000 compared to the y1 equation.

[Alireza Babaei]

Yes, I'm doing with 'format long', ...

here is my abstract: I will go agead with both cases (root1 = 0.14 and root1 = 1.87) to see what is the response of the system and then decide about which one to nominate as the rirst root

[Walter Roberson]

You are working with two different equations. You need to figure out which of the two equations is right.

I suspect that when you built

y = @(x)((10^-9)*x*(cosh(x)*sin(x)-sinh(x)*cos(x)) - ((10^-9)-1*(x^4)*(1+cosh(x)*cos(x)))) 

that you made a mistake.