Preben Alsholm

MaplePrimes Activity


These are replies submitted by Preben Alsholm

@janhardo That can be done:
 

DEtools:-DEplot(sys,[x(t),y(t)],t=0..20,[[x(0) = 1, y(0) = 0]],x=-1..1,y=-1..1,stepsize=0.1);
DEtools:-DEplot(sys,[x(t),y(t)],t=0..20,[[x(0) = 1, y(0) = 0]],x=-1..1,y=-1..1,stepsize=0.1,animatecurves=true);
DEtools:-DEplot(sys2,[x(t),y(t)],t=0..200,[[x(0) = 1, y(0) = 0]],x=-20..1,y=-6..6,stepsize=0.1);
DEtools:-DEplot(sys2,[x(t),y(t)],t=0..200,[[x(0) = 1, y(0) = 0]],x=-20..1,y=-6..6,stepsize=0.1,animatecurves=true);

Using sys and sys2 from my answer above.

@janhardo Here is another road to Rome for the simple example above:
 

restart;
eq1 := diff(x(t), t) = y(t);
eq2 := diff(y(t), t) = -x(t);
beginwaarden := x(0) = 1, y(0) = 0;
sol:=dsolve({eq1, eq2, beginwaarden}, {x(t), y(t)});

L:=subs(sol,[x(t),y(t)]); # List with value of x(t) first
XY:=unapply~(L,t);  # You can use any name on the left.
XY(tau);
plot(XY(t),t=0..2*Pi);
plot(XY,0..2*Pi);
## Now the alternative road:
L;
indets(L,function); # Just warming up.
XY2:=evalindets(L,function,s->op(0,s)); # The alternative
plot(XY2(t),t=0..2*Pi);
plot(XY2,0..2*Pi);

 

@janhardo Is this what you want to accomplish:
 

restart;
eq1 := diff(x(t), t) = y(t);
eq2 := diff(y(t), t) = -x(t);
beginwaarden := x(0) = 1, y(0) = 0;
sol:=dsolve({eq1, eq2, beginwaarden}, {x(t), y(t)});

L:=subs(sol,[x(t),y(t)]); # List with value of x(t) first
XY:=unapply~(L,t);  # You can use any name on the left.
XY(tau);
plot(XY(t),t=0..2*Pi);
plot(XY,0..2*Pi);

 

@janhardo Do this:

sol:=dsolve({eq1, eq2, beginwaarden}, {x(t), y(t)});
             
subs(sol,[x(t),y(t)]); # List with value of x(t) first

if you use rhs elementwise as in rhs~(sol) then the result is a set of the right hand sides of the elements in sol.
That leaves it to set ordering to determine the order, not you.

@nm I suspect that you already know this.
 

restart;
integrand:=(2*x^n+1)/(x^(n+1)+x);
res2021:=convert(eval(integrand,n=2021),parfrac): # Takes a while
length(res2021); #19646
simplify(res2021-eval(integrand,n=2021)); # 0

In fact it also works for n=2023..2030.
for n = 2030 it takes between 10 and 15 minutes.

.

@Thomas Richard I tried replacing the number 2022 by n:
 

restart;
integrand:=(2*x^2022+1)/(x^2023+x);
integrand_n:=(2*x^n+1)/(x^(n+1)+x);
infolevel[int]:=3:
int(integrand_n,x);
eval(%,n=2022);

That revealed that a Risch-Norman integrator was successful.
Continuing with:
 

int(integrand_n,x,method=_RETURNVERBOSE);

gave the information that norman, meijerg,and risch were successful.

@Kitonum Your workaround works in 2024.0 too and is very fast.

@ I  have Mathematica 13.2, but I use it extremely little.
Thus I cannot step in right away.

@Axel Vogt I noticed that you have an earlier build.
The results you get, however, are the same in RC4: `Maple 2024.0, X86 64 WINDOWS, Mar 01 2024, Build ID 1794891`

@Rouben Rostamian  That is very nice.

@ Why is the sequence B in inv_sierp being created?  It isn't used for anything.
B is a sequence of zeros and ones.

Another problem. K is a sequence of 0 and 2 only. So never odd.
Thus a part of the code can be commented out (as well as the B's):

for j to k do 
    km := K[k - j + 1]; # either 0 or 2
    t := t*2^d + cd - km; 
    (*if d = 2 then 
      if evalb(type(km, odd)) then t := 1 - t end if; # Never happens
    end if; 
    if evalf(t) < 0 then t := t + 1 end if; *)  # Never
  end do;

You get the same result as before after commenting this out.
### Note:

I tried using the original inv_sierp code with nothing commented out.
The only change to the procedure was that I made it return t ,[K], [B]
 

return t ,[K],[B]
end proc:
d := 2;
k := 20;
t := 0.6164;
P := sierp(d, k, t);

res:= inv_sierp(d, k, P);

res was:
res := 360798224641/2,
[2, 2, 2, 0, 2, 2, 2, 2, 2, 2, 0, 0, 2, 2, 2, 2, 0, 0, 0, 2],
[1, 0, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 1, 1, 1, 0, 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 1]

@ Yes, sierp now works beautifully.
I tried with N = 4, 5, 6, 7, and 8. It is very nice. For N=8 I changed the thickness to 1.

I'm now convinced that somehow this algorithm has to be iterated.
Try this:
 

restart;
with(Fractals:-LSystem);
with(LSystemExamples);
PlotExample(SierpinskiTriangle,7);
for i from 1 to 10 do PlotExample(SierpinskiTriangle,i) end do;

Notice that PlotExample(SierpinskiTriangle,1);  looks like this:

I tried

plot(P,view=[0..1,0..1]);

and got:

Comparing this with Fig. 1 and Fig. 2  in the pdf file obtained from your link made me wonder if we are misunderstanding the algorithm.

@ I tried putting in a return at this place in the code:

K := K, km;
  end do;
  return K;

What I got was:
 

N := 500; 
P := [seq(sierp(2, 20, j/N), j = 0 .. N)]:

numelems(P); ## 10020
{P[]}; # 6 members {-1, 0, 1, 2, 3, 5}
{P[1..1260][]};# {-1, 0}
{P[1261..3760][]};# {1}
{P[3761..6260][]};# {2}
P[1]
{P[6261..8760][]};# {3}
{P[8761..8762][]};  # {0}
{P[8763..10020][]};  # {0, 5}

Here we see 6 different results.
In your code you have this:

for j to k do
    km := K[k - j + 1];
etc

Since k in this case k is 20 only the first 20 of the 10020 K-values will be used:
Whether that is a problem I have no idea.

3 4 5 6 7 8 9 Last Page 5 of 229