## 19 Badges

8 years, 270 days

## Example...

A:=Matrix(20,2, (i,j)->`if`(j=1,i,(i-10)^2)):
B:=Matrix(20,2, (i,j)->`if`(j=1,i,0.75*(i-10)^2)):

PLOT(CURVES(A,LEGEND("1st"),COLOR(RGB,1,0,0)),
CURVES(B,LEGEND("2nd"),COLOR(RGB,0,1,0))    );

or

pA:=PLOT(CURVES(A,LEGEND("1st"),COLOR(RGB,1,0,0))):
pB:=PLOT(CURVES(B,LEGEND("2nd"),COLOR(RGB,0,1,0))):
plots:-display(pA,pB);

## singular...

Implicitplot algorithm needs sign changes in order to determine the points.
Mathematically, in this case, all the points of the surface (x+y+z)^2=0 are singular.

Of course, you can use the equivalent implicitplot3d(x+y+z=0, x=-5..5,y=-5..5,z=-5..5) to obtain the plane.

## Dancing with ... Maple...

Each of n dancers placed at random positions on the dance-floor randomly chooses one "friend" and one "enemy".
At each step, every dancer:
- moves 0.5% closer to the centre of the floor
- takes a large step towards his friend
- and a small step away from his enemy.
At random intervals the dancers re-choose their friends and enemies.

mu = number of steps between frames
rechoose = probability (%) to rechoose friends and enemies after mu*numframes steps.

n:=1000: numframes:=200: mu:=10: rechoose:=10:
x,y:='LinearAlgebra:-RandomVector(n,generator=-1.0 ..  1.0,datatype=float[8])' \$2:
FE:=proc() global F,E; #friend,enemy
F,E:='LinearAlgebra:-RandomVector(n,generator=1..n,datatype=integer[4])' \$2;
end:
FE():
step:=proc(x::Vector(datatype=float[8]),y::Vector(datatype=float[8]),
E::Vector(datatype=integer[4]),F::Vector(datatype=integer[4]))
option autocompile;
local i::integer[4],ex::float[8],ey::float[8],ed::float[8],
fx::float[8],fy::float[8],fd::float[8];
to mu do
for i to n do
ex:=x[E[i]]-x[i];ey:=y[E[i]]-y[i]; ed:=sqrt(ex^2+ey^2);
fx:=x[F[i]]-x[i];fy:=y[F[i]]-y[i]; fd:=sqrt(fx^2+fy^2);
x[i] := 0.995*x[i] - 0.01*ex/(ed+0.01) + 0.02*fx/(fd+0.01);
y[i] := 0.995*y[i] - 0.01*ey/(ed+0.01) + 0.02*fy/(fd+0.01);
od; od;
end:
for k to numframes do
p[k]:=plot(x,y);
step(x,y,E,F);
if rand(1..100)()<=rechoose then FE() fi;
od:
plots[display](seq(p[k],k=1..numframes),style=point,axes=none,insequence=true);

Edited. Added definition, comments and a few changes in code.

## incompatible...

Your system is incompatible. You have posted a similar system recently and you were told how to see this (e.g. remove the last equation).

## Iterator or cartprod...

allcomp:=proc(F, n)
local M:=Iterator:-CartesianProduct(F\$n);
[seq(`@`(entries(u,nolist)), u=M)];
end;

allcomp([sin, cos, t->t^2], 3)(x);  #example

For Maple<2016 use:
allcomp2015:=proc(F,n)
local T:=combinat:-cartprod([F\$n]),r:=NULL;
while not T[finished] do r:=r,`@`(op(T[nextvalue]()))  end do;
[r]
end;

## int...

To integrate f(x) over the interval [a,b] you can use

Int(f(x)*Heaviside(x-a)*Heaviside(b-x), x=-infinity..infinity);

But why would you want to do so, when Maple is able to find the antiderivative?

J:=int(c-sqrt(a+b*(v*x+u)^2), x);
x1,x2 := -(b*u+sqrt(b*c^2-a*b))/(b*v) , (-b*u+sqrt(b*c^2-a*b))/(b*v);
eval(J,x=x2) - eval(J,x=x1); simplify(%);

## Nothing to fix...

Your surface is actually a curve; it depends only on r = sqrt(u^2+v^2).
Near r=0 this curve is almost a line.

## int...

You want

int((((x/b)^c-d))*A*x^e, x=0..X) assuming c>0,e>0;

## Customized RRef...

@logiaco

For one parameter only, SmithForm e.g. could be used.
The only reasonable solution for >1 parameters I see
is to write a customized version of ReducedRowEchelonForm
e.g. RR(A, params:=indets(A))
having as result a sequence of lists

[cond, R]

where cond is  a list of conditions on parameters and R is the corresponding echelon form.

e.g.

[[Determinant(A)<>0], IdentityMatrix(n)], ...

This should not be difficult but I am not sure whether you really need it.

## Not compatible...

The system is not compatible. Execute in your worksheet:

solve({A[1], A[2], A[3], A[4], A[5]}, [a, b, c, d, e, f]);
[[a = 2+b, b = b, c = -2*b, d = -6+4*b, e = -4*b+3, f = -1+3*b]]
subs(%[], A[6]);
43 = 51

## OK for your corrected version...

eqx:=diff(f(b),b,b) = -(1/2)*(6*(diff(f(b), b))*f(b)-8*b*(diff(f(b), b))-2*f(b)+diff(f(b),b)+2)*(diff(f(b),b))/((-b+f(b))*(-1+f(b))):

ic:=f(3/8)=0, D(f)(3/8)=d:

s:=dsolve({eqx,ic},numeric,parameters=[d]):

pp:=proc(u)
s(parameters=[d=u]);
abs( eval(f(b),s(1/2)) - 1/2)
end:

N:=100000: V:=Vector(N/10, i-> pp(2.4+i/N)):

min(V):min[index](V):
s(parameters=[2.4+%/N]);
#  [d = 2.485270000]
s(3/8);
# [b = .375000000000000, f(b) = 0., diff(f(b), b) = 2.485270000]
s(1/2);
#[b = .500000000000000, f(b) = .499859411261823, diff(f(b), b) = 1.00693188008423]
plots[odeplot](s,3/8..1/2);

## using parameters...

The ODE may be expressed wrt f''(b) by:

eq:=diff(f(b), b, b) = (1/2)*((-8*b+6*f(b)+1)*(diff(f(b), b))^2+(-2*f(b)+2)*(diff(f(b), b))+2*b-2*f(b))/(f(b)*(b-f(b)));

Note that at both boundary values, the denominator is 0 (!).

Let's try using parameters for an IC problem:
b0:=(3/8+1/2)/2:
ss:=dsolve({eq, f(b0)=A, D(f)(b0)=B}, numeric, parameters=[A,B]):
pp:=proc(u,v)
ss(parameters=[A=u,B=v]);
abs(eval(f(b),ss(3/8)) - 0) + abs( eval(f(b),ss(1/2)) - 1/2)
end:

We want to minimize pp in order to be as close as possible to your conditions.
Here only a crude attempt is made:

M:=Matrix(37, (i,j) -> pp(i/100,j/100)):
min[index](M);min(M);

37, 37
0.44450458358379435

ss(parameters=[A=0.37,B=0.37]);
ss(3/8);
[b = .375000000000000, f(b) = .362588443481166, diff(f(b), b) = -0.179678246588966e-1]
ss(1/2);
[b = .500000000000000, f(b) = .418083859897371, diff(f(b), b) = 1.17994432455953]

So, it seems that it's not possible to satisfy acceptably your IC.

## Hidden definition...

Executing
eval(L);
in your document ==>
table([c1 = .5*L, c2 = .5*L])

So, you have defined somewhere a hidden infinite recursive definition for L.

That is why I always prefer a clean worksheet mode (with 1D math if possible)
such that everything is visible!

Edit.

Ok, looking closer, the recursive definition is not hidden. It appears at the beginning of the document:

L[c1] := 0.5*L;
L[c2] := 0.5*L;

So, just change:

L[c1] := 0.5*K;
L[c2] := 0.5*K;
and it works

## [0,1], [1,2]...

Maple computes correctly the integral in [0,1] and [1,2], but not in [0,2]. It's a pity.

J:=Int(floor(x^2), x = 0..2);
IntegrationTools:-Split(J,1): value(%);
5-sqrt(2)-sqrt(3)

## equivalent...

They are mathematically equivalent, but `simplify/D` does not implement this. But,

e1 := (D[1]@D[1])(f)(x);
convert(e1,diff): convert(%,D): lprint(%);
((D@@2)(f))(x)

 First 99 100 101 102 103 104 105 Last Page 101 of 113
﻿