Question: Where is the error in evalf or do?

This worksheet has a system of ode's ( a second order and a first order), a power series solution near 0 to use for inital conditions, and a function W of the vairables that I want to compute for certain values of x < π/2. W is a ratio of 2 positive values and that ratio seems to be negative when evalf is applied to the answer.!!! The true answer (without the evalf) are positive. Also, applying evalf to the correct answers does not produce those silly negative answers.

• >  restart;

• >  alias(B=b(x));alias(Bp=diff(b(x),x));alias(f=F(x));alias(fp=diff(F (x),x));alias(fpp=diff(F(x),x,x));

  B B, Bp

  B, Bp, f B, Bp, f, fp

  B, Bp, f, fp, fpp

• >  eq1 := sin(x)*(1+B*cos(x)^2)*fpp-(2*((f^2-2*B)*cos(x)^2-f^2+3*B* (1/2)-1))*cos(x)*fp+2*sin(x)*((f^2-3*B*(1/2))*cos(x)^2-f^2+(1/2)* B-1/2)*f = 0;eq2 := (-B*cos(x)^6+(B-1)*cos(x)^4+cos(x)^2)*fp^2+2*f* cos(x)*sin(x)*(cos(x)-1)*(cos(x)+1)*(1+B*cos(x)^2)*fp+cos(x)*sin(x) *Bp+B*cos(x)^6*f^2-(2*(B-1/2))*f^2*cos(x)^4+B*cos(x)^2*f^2-f^2+B = 0;

  eq1 sinx 1 Bcosx2fpp 2 f2 2Bcosx2 f2 3B 1cosxfp 2

  2sinx f2 3B cosx2 f2 B 1 f=0 222

  eq2 Bcos x 6 B 1 cos x 4 cos x 2 fp2 2fcos x sin x cos x
  1 cosx 1 1 Bcosx2 fp cosx sinx Bp Bcosx6f2 2 B

  1 2

(1)

f2cos x 4 Bcos x 2f2 f2 B=0 > f2p:=solve(eq1,fpp);

f2p

(2)

(3)

1 2sin x f3cos x 2 sin x 1 Bcos x 2

2cos x 3fpf2
3sinx fBcosx2 4cosx3fpB 2sinx f3 2cosx fpf2 sinx fB

3cosx fpB sinx f 2cosx fp

• >  fs:=x^8*(127/604800*C+653/17010*(C^5)+2/81*(C^7)+773/1020600*(C^3)) +x^6*(17/5670*(C^3)-4/81*(C^5)+31/15120*C)+x^4*((1/9)*C^3+7/360*C)+ (1/6)*x^2*C+C:

• >  Bs:=x^6*(20/189*(C^4)+4/315*(C^2))+4*x^4*C^2*(1/45)+2*x^2*C^2*(1/3) :

• >  Digits:=15:

• >  ic:=(c)->{F(1e-5)=subs([x=1e-5,C=c],fs),D(F)(1e-5)=subs([x=1e-5,C= c],diff(fs,x)),b(1e-5)=subs([x=1e-5,C=c],Bs)}:

> nans:=dsolve({eq1,eq2} union ic(c),{f,B},numeric,stiff=true,

(4)

> parameters=[c]);

nans proc x_rosenbrock ... end proc

(4)

(5) (6)

> W:=fp^2*cos(x)/f^2;

W fp2 cos x f2

> Wp:=simplify(subs(fpp=f2p,diff(W,x)));

cosxsinx 1 Bcosx2fp2 f27Bcosx4 f25B3cosx2fpcosxf2

1

1
Wp f3sin x 1 Bcos x 2 4fp 2 f 4

444

3B cosx2 f2 B 1 f2sinx 222

• >  k:='k': for k from 4 to 5 do nans(parameters=[k]); evalf(eval([W, Wp], nans(Pi/2-1e-10)),5); end do;

  c = 4.
  482.73, 1.3965 108

  c = 5.
  5.6838 106, 1.4948 1013

• >  k:='k': for k from 4 to 5 do nans(parameters=[k]); eval([W,Wp],

  nans(Pi/2-1e-10)); end do;

  c = 4. 0.0131412243081803, 1.31416436948280 108

  c = 5.
  154.729383663460, 1.54702786932958 1012

• >  evalf([0.0131412243081803, -1.31416436948280*10^8],5);

  0.013141, 1.3142 108

(7)

(8) (9)