## 12785 Reputation

8 years, 317 days

## periodic extension...

f:=x->piecewise(x<1,x^2,(2-x)^2):

f_per:=x->f(x-2*floor(x/2)):

## similar...

seq is faster but \$ works in symbolic contexts e.g.

[x[i], y[i]] \$ i=1..n;

diff( exp(a*x),  x\$n);

where n is a symbol.

## simplify symbolic...

You may use

Be aware that the validity depends on the domain of u in this case.

## 1D math...

Why don't you switch to 1D math and use simply:

diff(u, y, y, t);

## evaluation...

It is a Maple evaluation peculiarity. In sum, the argument  N[i] is evaluated first (before the index i is defined).

You must use

sum('N[i]', i = 1..4);
1

Or, better, use add, for which the evaluation rule is different.

1

## Change of variables (Rouben's idea)...

restart;
p:=10:
dx_:=seq(x[i]=0..1,i=1..p-1):
H:=int( Heaviside(xp-sx_)*Heaviside(1-xp+sx_),  dx_ ):
int( exp(-xp^3)*H, xp=0..p,  numeric);
0.000002546120390

# However the method is slow for higher accuracy and does not provide the exact result.
# It has the advantage that it works in other situations too.

## There are an infinity of equivalent syst...

There are an infinity of equivalent systems: e.g.

F(a, b, c) = F( s/RootOf(_Z^2-s^2+s), -RootOf(_Z^2-s^2+s)/s, RootOf(_Z^2-s^2+s) )

where F is any bijection R^3 --> R^3 (or C^3 --> C^3).

Here, an interesting situation is to eliminate the RootOf. For this, denote it by r, and use its definition ==>

[a = s/r, b = -r/s, c = r, r^2-s^2*r+s=0 ];
# Then use eliminate
eliminate(%,r);

So, your equivalent system without RootOfs is

## exact and numeric...

restart;
p:=10;
dx:= seq(x[i]=0..1,i=1..p);

We want to compute Int(exp(-h^3),  dx)
(in older versions of Maple, dx must be replaced by [dx]).

Taking the power series of exp, we need to compute Int( h^k,  dx).
Maple is able to compute it symbilically, but the direct computation is extremely slow (for p=10) and I had to interrupt and help Maple:

f:=u^k:

for j to p do
simplify(int(subs(u=u+t,%),t=0..1),size) * (k+j)
assuming k>=0,u>0 od: simplify(%):
Jp:=subs(u=0,%)/mul(k+j,j=1..p):
J:=unapply(Jp,k):
result:=sum(J(3*k)/k!*(-1)^k,k=0..infinity);

evalf[50](result);

0.000002546120389594240309918560289962734341705715230458

It was a pure luck that the series was obtained symbolically.
But note that anyway we can compute it numerically.
It is interesting to note that even if the convergence is very fast, the first terms are huge.
e.g.
max(seq(J(3*k)/k!,k=0..3000)):evalf(%);

max(seq(J(3*k)/k!,k=3000..4000)):evalf(%);

## exact...

The integral can be computed exactly: just replace Int by int.

For a 6 fold integral, a standard numerical integration is almost impossible. Consider the fact that if one has 100 division points in each interval, the function will be evaluated 10^12 times!

## series...

f:=series( exp(k*t)*cos(w*t),t=0);
series(f/exp(k*t),t=0);

## Eigenvectors...

For the 2nd sheet you should use Eigenvectors which is designed exactly for this purpose.

A:=eval(q,omega=0);
C:=eval(A-q,omega^2=1);
Eigenvectors(A, C);

(the eigenvalues correspond now for lambda = omega^2)

## df/dv...

If p is not available, then use:

f:=3*(r/sqrt(a)) + (r/sqrt(a))^2:
v:= r/sqrt(a):
diff(f,r)/diff(v,r);

## Y(arguments)...

You must change Y with Y(r1,r2,r3) and it will work.
Also, the inverse transforms are not necessary (being very simple), so it's enough:

PDEtools[dchange]({r1=v/2+w/4, r2=u/2+w/4, r3=u/2+v/2},  del_1 + del_2);

## solution...

The system is very special. If we try solve, we get:

sys := {-6-4*y-x-(1+y)*x+sqrt((4*(1+y))*(2+x)*(4+2*y+x)+(-(1+y)*x+2+x)^2), (2*(4+2*y+x))*(1+y)-(1+y)*x+2+x+sqrt((4*(1+y))*(2+x)*(4+2*y+x)+(-(1+y)*x+2+x)^2)-(2+y)*(-(1+y)*x+2+x+sqrt((4*(1+y))*(2+x)*(4+2*y+x)+(-(1+y)*x+2+x)^2))};

solve(sys);

This means that the system reduces to an identity in a certain region of the plane (x,y); a warning would have been very useful!.
Unfortunately, in such cases:
- Maple in unable to find this region
- fsolve does not work.

So, we will have to maniputate the system, or help Maple.

simplify(sys) assuming x*y+2*x+4*y+6>=0;
{0}
simplify(sys) assuming x*y+2*x+4*y+6<0;

solve(sys) assuming x*y+2*x+4*y+6>=0;

solve(sys) assuming x*y+2*x+4*y+6<0;

Conclusion: the system reduces to identities in the region x*y+2*x+4*y+6>=0;
(it is the "interior" of a hyperbola):

In the "exterior" of this hyperbola, the solution is on a curve (also a hyperbola). You may want to plot it.

## Fourier...

For the function x : [0,1] --> R you probably want the Fourier expansion on each of the N intervals [(n-1)/N, n/N], 1<=n<=N.

In this case, X(m,n) must be the Fourier coefficient associated to the restriction of x() to this interval.

 First 106 107 108 109 110 111 112 Page 108 of 114
﻿