how to solve for non linear system of equations containing series terms.

Hey all, I want to solve for following non linear system of equations:

I have only 3 unknown: k=x(1); F=x(2); epsilon= x(3); And the corresponding nonlinear system of equations is :

Equation 1: 0.94=[(k^2+1)*epsilon-2*F]/[(k^2-1)*epsilon]
Equation 2: F=(pi/2)*summation n varies from 0 to inf {[(2n!)/(2^(2n)*(n!)^2)]^2}*[k^(2n)]
Equation 3: epsilon=(pi/2)*summation n varies from 0 to inf {[(2n!)/(2^(2n)*(n!)^2)]^2}*[k^(2n)]/(1-2n)

To solve for this I have written following code:

function F = myfun(x)
curvediff=0.94;
n=0:1:10e3;
N1=factorial(2.*n);
D1=2.^(2.*n);
D2=factorial(n).^2;
ND=N1./(D1.*D2);
PI=ND.*ND;
PII=PI./(1-2.*n);
EI=PI.*x(1).^(2.*n);
EII=PII.*x(1).^(2.*n);
F=[curvediff-((x(1)*x(1)+1)*x(3)-2*x(2))/((x(1)*x(1)-1)*x(3));x(2)-pi/2.*sum(EI);
  x(3)-pi/2.*sum(EII)];

And when I try to solve for this I am getting following error :

>> x0=[9.347;5;1.02];
>> x=fsolve(@myfun,x0);
??? Error using ==> trustnleqn at 28
Objective function is returning undefined values at
initial point. FSOLVE cannot continue.

Error in ==> fsolve at 366
  [x,FVAL,JACOB,EXITFLAG,OUTPUT,msgData]=...

Kindly help me to debug this error. Thanks in advance.

Nikhil

Answer 1

[Roger Stafford]

Nikhil, your method of evaluating the two infinite series will not work for the sizes of the values in 'n'. Factorial(2e4) is an indescribably large value far, far outside the range of matlab's double precision numbers. You should be evaluating the quantities in your series iteratively in terms of their previous values in 'n'. For example the term in the series for F with n = 3 is:

((1*3*5)/(2*4*6))^2*k^6

To get to the next term with n = 4, the above can be multiplied by

 (7/8)^2*k^2

to arrive at

 ((1*3*5*7)/(2*4*6*8))^2*k^8

Generalizing this technique you can avoid the horrendous sizes of the factorials when n becomes large.

(Note: I am assuming you meant {[(2n)!/(2^(2n)*(n!)^2)]^2}*[k^(2n)] in your expression, as in your previous posting.)

Answer 2

[Roger Stafford]

Based on the same principle as in my other answer you can compute the sum of those series using 'sum' and 'cumprod' in a manner which also avoids computing factorials. For example

 n=0:1:(some large integer);
 n1 = n(2:end);
 F = pi/2*sum(cumprod([1,((2*n1-1)./(2*n1)*k).^2]));

(The 'F' here refers to the quantity in Equation 2, not the F in 'myfun'.)

(Note that the quantity 'k' must be kept less than 1. Otherwise the series diverges to infinity.)