How to set eps^n =0 for all n>= 10 where eps is a symbolic variable.

Here is some Reduce code for an perturbation calculation:

1 factor eps;

2 x:=-1;

3 let eps^10=>0;

4 for iter:=1:2 do begin

5 write res:=x^2+x-eps;

6 write x:=x+a*res;

7 end;

and corresponding Mathematica code:

1 Clear[x,res,eps]

2 eps^n_ ^:=0 /; n>6;

3 x=-1;

4 While[res=!=0,

5 res=x^2+x-eps // Simplify;

6 x=Expand[x+res];

7 Print[x]

8 ];

I would like to know how to do line 3 in the Reduce code that sets all powers of eps^10 or greater equal to zero. The Mathematica code on line 2 does the same. I thought I could get subsex to do this but I can not figure out how to make say eps^12 =0 if eps^10=0, ie. obey the algebraic rules for powers like reduce does: eps^12 = eps^10 eps^2 = 0 since we let eps^10=0. I think maple's function algsubs might be able to do this.

[Walter Roberson]

It is not clear that you are trying to do this work in MATLAB. subsex does exist inside of MuPAD but it is very rare that you would use evalin(symengine) or feval(symengine) to deliberately call a MuPAD operator.

[John D'Errico]

It would have been a lot easier if you at least started to write MATLAB code. Instead, we need to figure out what you are doing from other languages, on a MATLAB forum. I may have gotten this right, but probably not.

Anyway, never use eps as a variable in MATLAB. eps is far too useful already, so I'll call it EPS. I'll guess this is what you want to do. I'll just run the iterations out a little ways to see what is happening.

syms EPS
x=-1;
for iter = 1:6
  res = expand(x^2+x-EPS)
  x = simplify(x+res);
end

res = 

res = 

res = 

res = 

res = 

res = 

So we have a simple polynomial in EPS at each iteration in res. After 4 iterations, there were no too large powers.

As well, I see the lowest power of EPS is always iter. So at each iteration, that lowest power of EPS increases by 1.

Next, we should recognize that if EPS^10 or greater is zero, then we can just leave those terms in the expansion. zero is zero, and adding or multiplying by some power of EPS that is zero does nothing more than generate higher powers of zero.

Or, we could put a line inside the loop, where we use coeffs to extract the individual terms of this polynomial. Then zero out the coefficients of all terms with an exponent at least 10, and recombine to get the desired truncated polynomial.

Regardless, everything will die at iteration 10.

for iter = 7:10
  res = expand(x^2+x-EPS)
  x = simplify(x+res);
end

res = 

res = 

res = 

res = 

[C,T] = coeffs(res)

C = 

T = 

Anyway, it is not clear where you are going with this, but it should not be even at all difficult. As I said, at each step, if you really want to kill off those higher powers of EPS, you can do it. Insider the loop, just zero out those coefficients, then recombine what is left in C and T into a new version, one that has no higher order powers of EPS.

[Jack]

I was not very clear as this is my first post. There is nothing special about eps, here is another/different example of generating a power series by Picard type iterations

syms x
y=1-2*x
for n=1:3
    y=1-2*x+6*int(int(y^2,x),x)
    n
end

I would like to have only powers of x^7 and lower so somehow in particluar a series that is O(x^8). if

Here is a reduce code that does this

1 let x^8=>0;

2 y:=1-2*x;

3 for n:=1:5 do write

4 y:=1-2*x+6*int(int(y^2,x),x);

And the output is for the first 5 iterations:

The reduce code line 1 makes the power series generated to be O(x^8). This is what I cannot figure out how to do in Matlab.

Answer 1

[Walter Roberson]

mapSymType(x, 'power', @(X) piecewise(children(X,1)==eps & children(X,2)>10, 0, X))

However, eps is normally a function in MATLAB that returns epsilon of its given parameter. eps by itself is normally equivalent to eps() which in turn is equivalent to eps(1) -- the epsilon of 1.0 . Normally the function is immediately evaluated, and eps to a non-symbolic power would normally be evaluated in double precision, rather than being held as a eps^n form. You would need to do something like

eps = sym('eps');

to get eps to be a symbolic variable that eps^10 would be represented symbolically.

[Jack]

I was not very clear as this is my first post. There is nothing special about eps, here is another/different example of generating a power series by Picard type iterations

syms x
y=1-2*x

y = 

for n=1:3
    y=1-2*x+6*int(int(y^2,x),x)
    n
end

y = 

n = 1

y = 

n = 2

y = 

n = 3

I would like to have only powers of x^7 and lower so somehow in particluar a series that is O(x^8). if

Here is a reduce code that does this

1 let x^8=>0;

2 y:=1-2*x;

3 for n:=1:5 do write

4 y:=1-2*x+6*int(int(y^2,x),x);

And the output is for the first 5 iterations:

The reduce code line 1 makes the power series generated to be O(x^8). This is what I cannot figure out how to do in Matlab.

[Walter Roberson]

syms x
y=1-2*x

y = 

for n=1:3
    y=1-2*x+6*int(int(y^2,x),x)
end

y = 

y = 

y = 

mapSymType(y, 'power', @(X) piecewise(children(X,1)==x & children(X,2)>=10, 0, X))

ans = 

[Jack]

Looks like this works, Thanks!

syms x
y=1-2*x

y = 

for n=1:6 
    y=1-2*x+6*int(int(y^2,x),x);
    y=mapSymType(y, 'power', @(X) piecewise(children(X,1)==x & children(X,2)>=8, 0, X))
    n
end

y = 

n = 1

y = 

n = 2

y = 

n = 3

y = 

n = 4

y = 

n = 5

y = 

n = 6

[Paul]

@Jack Dockery

Keep in mind that int only returns one particular antiderivative, not a general expression for the antiderivative.

For example

syms x
y=1-2*x

y = 

After the first iteration

y = 1-2*x+6*int(int(y^2,x),x)

y = 

But if we add the constants of integration explicitly, we'd have

y = 1 - 2*x

y = 

syms C1 C2
y = 1-2*x+6*int(int(y^2,x) + C1,x) + C2

y = 

Also, sometimes int can return an unexpected form of an antiderivative depending on the integrand.

Just some things to keep in mind if you care about the particular form of the antiderivative as the code proceeds through the iterations.

@Walter Roberson

Does mapSymType automatically expand the first input argument? If not, maybe it would miss something? For example, after the first iteration on n the expression y is a function of x^4, but would mapSymType see that if the goal was to eliminate powers greater than 3?

Also, maybe it would be better to do the mapSymType once after the iteration IF the iterated expressions were of a form where a high positive power of x could multiply a large negative power of x resulting in a term in x with a low power that's desired to keep. I guess that depends on the original expression for y and the iteration expression.

[Walter Roberson]

@Paul

mapSymType works with the internal form of the expression as-is. It does not expand() the expression at all.

It might be best to add an expand()

y = mapSymType(expand(y), 'power', @(X) piecewise(children(X,1)==x & children(X,2)>=8, 0, X))

[Jack]

Thanks for the help. The constants of integration have been choosen so that the solution to the ODE

y''=6 y^2

satisfies

y(0)=1 and y'(0)=-2

this is why we have the 1-2x plus the iterated integrals. Again sorry for the confusion.

[Paul]

No confusion.

Keep in mind the potential problem going the other way, i.e., int is allowed to return a result with an additive constant, which would affect the iterates. I don't know that it will for any case you're interested, but it's allowed, so that might be something to check for and deal with as needed after each iteration.

Just out curiosity ... I know next to nothing about Picard iteration, but I thought that approach is based on solving a definite integral at each iteration, not an antiderivative.

[Walter Roberson]

int is allowed to return a result with an additive constant

It is not especially common that int() returns results with additive constant, but it happens from time to time. In particular, int() detects some trig identities and rewrites them and then integrates the result, so you might get a result that has an implied "plus 1" in it (for example)

I do not have an example at the moment.