The goal here is to determine how many roots there are in a nonlinear function we are interested in solving. For this example, we use a cubic polynomial because we know there are three roots.
clc; clear all; close all;
This ony works for a polynomial, it does not work for any other nonlinear function.
roots([1 6 -4 -24])
ans = -6.0000 2.0000 -2.0000
x = linspace(-8,4); fx = x.^3 + 6*x.^2 - 4*x - 24; plot(x,fx)
count the number of times the sign changes in the interval. We determine the sign of each element (-1, 0 or +1) and then use a a loop to go through each element. A sign change is detected when an element is not equal to the previous element. Note this is not foolproof; we may over count roots if we get a sequence -1 0 1. We check for that by making sure neither of the two elements being compared are zero.
sign_fx = sign(fx); count = 0; for i=2:length(sign_fx) if sign_fx(i) ~= sign_fx(i-1) && sign_fx(i) ~= 0 && sign_fx(i-1) ~= 0 count = count + 1; disp(sprintf('A root is between x = %f and %f',x(i-1),x(i))) end end disp(sprintf('The number of sign changes is %i.',count))
A root is between x = -6.060606 and -5.939394 A root is between x = -2.060606 and -1.939394 A root is between x = 1.939394 and 2.060606 The number of sign changes is 3.
this one line code also counts sign changes. the diff command makes all adjacent, identical elements equal zero, and sign changes = 2 (1 - (-1)) or (-1 - 1). Then we just sum up the diff vector, divide by two, and that is the number of roots. Note that no check for 0 is needed; you will get 1/2 a sign change from -1 to 0. and another half from 0 to 1. They add up to a whole sign change.
ans = 3
Method 1 is functional but requires some programming. Method 1a is concise but probably tricky to remember.
using events in an ODE solver Matlab can identify events in the solution to an ODE, for example, when a function has a certain value, e.g. f(x) = 0. We can take advantage of this to find the roots and number of roots in this case.
with f(-8) = -120
xspan = [-8 4]; fp0 = -120; options = odeset('Events',@event_function); [X,Y,XE,VE,IE] = ode45(@myode,xspan,fp0,options); XE % these are the x-values where the function equals zero nroots = length(XE); disp(sprintf('Found %i roots in the x interval',nroots)) % lets compare the solution to the original function to make sure they are % the same. figure plot(X,Y,x,fx) legend 'numerical' 'analytical'
function dfdx = myode(x,f) dfdx = 3*x^2 + 12*x - 4; function [value,isterminal,direction] = event_function(t,y) value = y - 0; % then value = 0, an event is triggered isterminal = 0; % do not terminate after the first event direction = 0; % get all the zeros % categories: ODEs, nonlinear algebra
XE = -6.0000 -2.0000 2.0000 Found 3 roots in the x interval