## 13613 Reputation

19 years, 253 days

## Undocumented?...

Algebraic:-RecursiveDensePolynomials seems to be undocumented. No help available?

## Undocumented?...

Algebraic:-RecursiveDensePolynomials seems to be undocumented. No help available?

## type name...

@PatrickT My reason for the type check X::name was that the procedure is not defined to handle vector input as argument number 2. Thus you should be told (in the form of an error message) if you try to pass any other object than a name.

The reason for the line

S:=select(has,indets(F(X),indexed),X);

was to allow functions containing indexed names other than the variables, e.g.

H:=X->Vector([a[8]*f(X[1],X[2]),g(X[1],X[2]),h(X[1],X[2])]);

where a[8] is just a constant.

## type name...

@PatrickT My reason for the type check X::name was that the procedure is not defined to handle vector input as argument number 2. Thus you should be told (in the form of an error message) if you try to pass any other object than a name.

The reason for the line

S:=select(has,indets(F(X),indexed),X);

was to allow functions containing indexed names other than the variables, e.g.

H:=X->Vector([a[8]*f(X[1],X[2]),g(X[1],X[2]),h(X[1],X[2])]);

where a[8] is just a constant.

## Problems with ode3 and the initial condi...

@goli (1) You forgot to include ode3 in the dsolve command.

(2) The right hand side of ode3 is undefined at z = 0.

(3) You require A(0) = 0, but A(z) appears in the denominator in ode3.

You can solve for H,L,R as before. Dont include A(z). The equation for A(z) you solve exactly (i. e. not numerically).

restart;
with(plots):
eq := z-> (H^2+(-1)*.27*(1+z)^3-(1/20000)*(1+z)^4)/(H^2)^.1 = .7299500000;
yp := implicitdiff(eq(z), H, z);

ode:=diff(H(z),z)=subs(H=H(z),yp);

ode1:=diff(L(z),z)=L(z)/(1+z)+(1+z)/H(z);

ode2:=diff(R(z),z)=1/H(z);

ode3:=diff(A(z),z)=(-2*A(z))/(3*z)+(beta^(1/3))/(3*z*H(z)*(A(z)^(1/2)));

E:=dsolve({ode3,A(0)=0});
p := dsolve({ode,ode1,ode2,L(0)=0,H(0)=1,R(0)=0}, numeric, output=listprocedure);

Hp,Lp,Rp:=op(subs(p,[H(z),L(z),R(z)]));
map(expand,rhs(E));
evalindets(%,specfunc(anything,Int),x->Rp(z));
AA:=eval(%,beta=Hp(0.35));
plot(AA,z=0..10);

## Problems with ode3 and the initial condi...

@goli (1) You forgot to include ode3 in the dsolve command.

(2) The right hand side of ode3 is undefined at z = 0.

(3) You require A(0) = 0, but A(z) appears in the denominator in ode3.

You can solve for H,L,R as before. Dont include A(z). The equation for A(z) you solve exactly (i. e. not numerically).

restart;
with(plots):
eq := z-> (H^2+(-1)*.27*(1+z)^3-(1/20000)*(1+z)^4)/(H^2)^.1 = .7299500000;
yp := implicitdiff(eq(z), H, z);

ode:=diff(H(z),z)=subs(H=H(z),yp);

ode1:=diff(L(z),z)=L(z)/(1+z)+(1+z)/H(z);

ode2:=diff(R(z),z)=1/H(z);

ode3:=diff(A(z),z)=(-2*A(z))/(3*z)+(beta^(1/3))/(3*z*H(z)*(A(z)^(1/2)));

E:=dsolve({ode3,A(0)=0});
p := dsolve({ode,ode1,ode2,L(0)=0,H(0)=1,R(0)=0}, numeric, output=listprocedure);

Hp,Lp,Rp:=op(subs(p,[H(z),L(z),R(z)]));
map(expand,rhs(E));
evalindets(%,specfunc(anything,Int),x->Rp(z));
AA:=eval(%,beta=Hp(0.35));
plot(AA,z=0..10);

## Yet another way...

@PatrickT

The following version has the advantage that it also accepts matrices of arbitrary dimensions mxn.

TaylorExpand := proc(F,X::name,n::posint) local d,k,X1;
d:=[op(1,F(X))][1];
X1:=op(0,F(X))(op(1,F(X)),subsop~(2=[seq(X[k],k=1..d)],op(2,F(X))),op(3..-1,F(X)));
mtaylor~(F(X),X1,n)
end proc:
TaylorExpand(M,X,2);
M2 := X-> Matrix(2,3,(i,j)->f[i,j](X[1],X[2]));
M2(X);
TaylorExpand(M2,X,2);
TaylorExpand(G,X,4);

And allowing for the number of variables to be different from the row dimension:

TaylorExpand := proc(F,X::name,n::posint) local S,nv,d,k,X1;
S:=select(has,indets(F(X),indexed),X);
nv:=max(op(op~(S)));
d:=[op(1,F(X))][1];
X1:=op(0,F(X))(op(1,F(X)),subsop~(2=[seq(X[k],k=1..nv)],op(2,F(X))),op(3..-1,F(X)));
mtaylor~(F(X),X1,n)
end proc:

H:=X->Vector([f(X[1],X[2]),g(X[1],X[2]),h(X[1],X[2])]);
TaylorExpand(H,X,4);

## Yet another way...

@PatrickT

The following version has the advantage that it also accepts matrices of arbitrary dimensions mxn.

TaylorExpand := proc(F,X::name,n::posint) local d,k,X1;
d:=[op(1,F(X))][1];
X1:=op(0,F(X))(op(1,F(X)),subsop~(2=[seq(X[k],k=1..d)],op(2,F(X))),op(3..-1,F(X)));
mtaylor~(F(X),X1,n)
end proc:
TaylorExpand(M,X,2);
M2 := X-> Matrix(2,3,(i,j)->f[i,j](X[1],X[2]));
M2(X);
TaylorExpand(M2,X,2);
TaylorExpand(G,X,4);

And allowing for the number of variables to be different from the row dimension:

TaylorExpand := proc(F,X::name,n::posint) local S,nv,d,k,X1;
S:=select(has,indets(F(X),indexed),X);
nv:=max(op(op~(S)));
d:=[op(1,F(X))][1];
X1:=op(0,F(X))(op(1,F(X)),subsop~(2=[seq(X[k],k=1..nv)],op(2,F(X))),op(3..-1,F(X)));
mtaylor~(F(X),X1,n)
end proc:

H:=X->Vector([f(X[1],X[2]),g(X[1],X[2]),h(X[1],X[2])]);
TaylorExpand(H,X,4);

## One way...

@PatrickT You could do like this, where the previous definition of d has been replaced by d:=[op(1,F(X))][1];

A few other changes have been necessary.

TaylorExpand := proc(F,X::name,n::posint) local d,q,i,j,X1,X2;
d:=[op(1,F(X))][1];
if type(F(X),Matrix) and not op(1,F(X))[2]=1 then error "Matrix must have one column only" end if;
X1:=Vector(d,i->[seq(X[k],k=1..d)]);
X2:=Matrix(d,1,(i,j)->[seq(X[k],k=1..d)]);
if not type(F(X),Matrix) then mtaylor~(F(X),X1,n) else mtaylor~(F(X),X2,n) end if
end proc:
M := X-> Matrix(2,1,(i,j)->f[i,j](X[1],X[2]));

TaylorExpand(M,X,3);

I understand that an assumption is that d is also the number of variables?

## One way...

@PatrickT You could do like this, where the previous definition of d has been replaced by d:=[op(1,F(X))][1];

A few other changes have been necessary.

TaylorExpand := proc(F,X::name,n::posint) local d,q,i,j,X1,X2;
d:=[op(1,F(X))][1];
if type(F(X),Matrix) and not op(1,F(X))[2]=1 then error "Matrix must have one column only" end if;
X1:=Vector(d,i->[seq(X[k],k=1..d)]);
X2:=Matrix(d,1,(i,j)->[seq(X[k],k=1..d)]);
if not type(F(X),Matrix) then mtaylor~(F(X),X1,n) else mtaylor~(F(X),X2,n) end if
end proc:
M := X-> Matrix(2,1,(i,j)->f[i,j](X[1],X[2]));

TaylorExpand(M,X,3);

I understand that an assumption is that d is also the number of variables?

## Elementwise and op...

@PatrickT Using one of the ways of finding the dimension of a vector given by acer you could do like this

TaylorExpand := proc(F,X::name,n::posint) local d,i;
d:=op(1,F(X));
mtaylor~(F(X),[[seq(X[i],i=1..d)]\$d],n)
end proc;

## Elementwise and op...

@PatrickT Using one of the ways of finding the dimension of a vector given by acer you could do like this

TaylorExpand := proc(F,X::name,n::posint) local d,i;
d:=op(1,F(X));
mtaylor~(F(X),[[seq(X[i],i=1..d)]\$d],n)
end proc;

## Ineffective?...

As much as I hate 2d-input I appreciate the effort to make input look like mathematical expressions the way they are usually written.

Thus I appreciate very much that you can define a function f, so that f(2) returns the value of f at 2. In your comment are you somehow trying to say that the use of procedures is ineffective or slow and that the use of eval is faster or more effective?

(I should note that I too have seen students define procedures in situations where they never make use of any other input than, say x.)

## Ineffective?...

As much as I hate 2d-input I appreciate the effort to make input look like mathematical expressions the way they are usually written.

Thus I appreciate very much that you can define a function f, so that f(2) returns the value of f at 2. In your comment are you somehow trying to say that the use of procedures is ineffective or slow and that the use of eval is faster or more effective?

(I should note that I too have seen students define procedures in situations where they never make use of any other input than, say x.)

## Homework problem?...

Your formulation suggests that this is a homework problem or similar. Am I mistaken?

 First 214 215 216 217 218 219 220 Last Page 216 of 229
﻿