Finding the root of a function

5 views (last 30 days)
Daenerys
Daenerys on 11 Dec 2017
Edited: Daenerys on 10 Apr 2018
I'm trying to find x in the following equation using the fzero function.
I'm not sure whether I should use int or integral, and my code doesn't work. Any help?
My code:
N = 100;
S1 = 0;
S2 = 0;
for n = 0:N
S1 = S1 + ((-1)^n/((n + 1/2)*pi)^4*tanh((n + 1/2)*pi*(1+x+x.^2)/2));
S2 = S2 + (1/((n + 1/2)*pi)^5*tanh((n + 1/2)*pi*(1+x+x.^2)/2));
end
cs = 1/2 - 4./(1+x+x.^2).*S1;
cp = 1/3 - 4./(1+x+x.^2).*S2;
syms x
fun = @(x) cs./cp;
q = int(fun,0,x)
gx = @(x) q - 0.0062;
x0 = [0 10];
x = fzero(gx,x0);
The error:
Error using subsindex
Function 'subsindex' is not defined for values of class 'function_handle'.
Error in line 69
q = int(fun,0,x)
  2 Comments
Geoff Hayes
Geoff Hayes on 11 Dec 2017
Lilach - please clarify what you mean by your code is not working. Is there an error? If so, please copy and paste the full error message here. Or, are you not getting the expected answer?
Daenerys
Daenerys on 11 Dec 2017
Thanks for the input, I've added my error.

Sign in to comment.

Accepted Answer

David Goodmanson
David Goodmanson on 12 Dec 2017
Edited: David Goodmanson on 12 Dec 2017
Hi Lilach,
After moving the syms x statement to the top so that the code runs, it is not so clear that it is going to get there. Here is a slightly different method.
A = .0062;
N = 100;
% integrate in z from 0 to x, subtract A, find root
p = fzero(@(x) integral(@(z) ratfun(z,N),0,x)-A, [0 .1])
% take a look at the integrand
z = 0:.001:1;
plot(z,ratfun(z,N))
function y = ratfun(x,N)
S1 = 0;
S2 = 0;
for n = 0:N
S1 = S1 + ((-1)^n/((n + 1/2)*pi)^4*tanh((n + 1/2)*pi*(1+x+x.^2)/2));
S2 = S2 + ( 1/((n + 1/2)*pi)^5*tanh((n + 1/2)*pi*(1+x+x.^2)/2));
end
cs = 1/2 - (4./(1+x+x.^2)).*S1;
cp = 1/3 - (4./(1+x+x.^2)).*S2;
y = cs./cp;
end
The result is p = 0.004647. From the plot, the integrand starts out at about 1.34, and the value of the integral is so small, .0062, that you would not expect the integrand to change much from 1.34. So you would expect the rectangle result p = .0062/1.34 which is pretty close.
  1 Comment
Daenerys
Daenerys on 12 Dec 2017
Thanks David, that method works perfectly!

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!