## 11749 Reputation

15 years, 295 days

## Maybe this...

This unassigns all the numbered members:

```restart;
with(Statistics):
X := a -> RandomVariable(Exponential(a)):
##
for n from 1 to 100 do
S := X(n);
# here are some operations to do
end do:
##
A:={anames(alluser)}:
IA:=indets(A,suffixed(_ProbabilityDistribution,integer)):
##
eval(_ProbabilityDistribution5);
exports(_ProbabilityDistribution5);
##
unassign(op(IA));
eval(_ProbabilityDistribution5);
```

But you might want to unassign all but the last one:

```N:=100:
for n from 1 to N do
S := X(n);
# here are some operations to do
if n = N-1 then
A:={anames(alluser)};
IA:=indets(A,suffixed(_ProbabilityDistribution,integer));
unassign(op(IA));
end if
end do:
anames(alluser);
```

## Infinitely many solutions...

Your integral equation has infinitely many solutions because the corresponding homogeneous equation has a one parameter family of solutions.
A particular solution can be found as vv did.

You will find my reasoning in the attached worksheet:

Preben Alsholm, March 17 2020.

 > restart;

Notice that I have replaced the coefficient 1/2 of the integral with the parameter a:

 > eq := Z(x) =  3/2-9/2*exp(-2*x)-9/2*exp(-x)+a*int(exp(-y)*Z(x-y), y= 0..ln(2));
 > EQ:=IntegrationTools:-Change(eq, x-y=t, t);

The corresponding homogeneous equation:

 > EQ0:=Z(x) = a*Int(exp(t - x)*Z(t), t = x - ln(2) .. x);

We find a particular solution to the inhomogeneous equation following vv:

 > ZZ:=x -> A*exp(-2*x)+B*exp(-x)+C;
 > expand(eval( (lhs-rhs)(EQ), Z=ZZ));
 > coeffs(%, [exp(x),exp(2*x)]);
 > ABC:=solve([%],{A,B,C});
 > Zsol:=eval(ZZ(x),ABC);
 > simplify(eval( (lhs-rhs)(EQ), Z=unapply(Zsol,x))); # Check

Now find solutions of the homogeneous equation by first differetiating it and thereby getting a delay differential equation:

 > INT:=combine(-1*(rhs=lhs)(EQ0)); # For use below
 > dde01:=(combine@expand)(diff(EQ0,x));

The delay differential equation:

 > dde0:=eval(dde01,INT);

The equation dde0 has as one solution Z(x) = exp(-x) for all values of a:

 > eval((lhs-rhs)(dde0),Z=(x->exp(-x)));
 > simplify(%);

but notice that Z(x) = exp(-x)  doesn't satisfy  the homogeneous integral equation unless a = 1/ln(2):

 > eval((lhs-rhs)(EQ0),Z=(x->exp(-x)));
 > factor(value(%));

Trying to find other solutions to dde0 of the form Z(x) = exp(r*x):

 > value(eval((lhs-rhs)(dde0),Z=(x->exp(r*x))));
 > u:=(combine@expand)(%/exp(r*x)); # Needs to be zero
 > eval(u,r=-1);

Numerically we see here that for a < 0 there is one solution for r and for a > 0 two solutions except at (as it turns out) at a = 1/ln(2).

One of the solutions is -1.

 > plots:-animate(plot,[u,r=-8..2,-1..1,labels=[r,"u"]],a=-2..3,frames=100);
 > usol:=expand(solve(u=0,r,AllSolutions));
 > zz:=op(indets(usol,suffixed(_Z,integer)));

We take the two real branches:

 > rm1:=eval(usol,zz=-1);
 > r0:=eval(usol,zz=0);
 > plot([rm1,r0],a=-2..3,color=[red,"Green"],thickness=3,legend=["rm1","r0"]);

The two branches meet at a = 1/ln(2):

 > eval(rm1-r0,a=1/ln(2));

Checking that the solutions Z(x) = exp(rrm1*x) and  Z(x) = exp(r0*x) actually satisfy dde0:

First rm1:

 > rr:=expand(eval((lhs-rhs)(value(dde0)),Z=(x->exp(rm1*x))));
 > simplify(rr);
 > rr:=expand(eval((lhs-rhs)(value(dde0)),Z=(x->exp(r0*x))));
 > simplify(rr) ;

Now checking that they also satisfy EQ0, but remember that actually exp(-x) is not a solution unless a = 1/ln(2):

 > rr:=expand(eval((lhs-rhs)(value(EQ0)),Z=unapply(exp(rm1*x),x)));
 > simplify(rr); # I see!  It should only be zero for 0 < a <= 1/ln(2)

Testing with a = 1/2:

 > evalf(eval(rm1,a=1/2));
 > rrc:=expand(eval[recurse]((lhs-rhs)(value(EQ0)),{a=1/2,Z=unapply(exp(rm1*x),x)}));
 > simplify(%); # 0
 > evalf(rrc);  # 0.
 > rrcn:=normal(rrc);
 > simplify(denom(rrcn));
 > evalf(%); # -0.825...
 > simplify(numer(rrcn));

We look at r0.

 > rr:=expand(eval((lhs-rhs)(value(EQ0)),Z=unapply(exp(r0*x),x)));
 > simplify(rr); # # It should only be zero for 1/ln(2) <= a

We look at a = 1/2:

 > eval(r0,a=1/2);
 > evalf(%);
 > rrc:=expand(eval[recurse]((lhs-rhs)(value(EQ0)),{a=1/2,Z=unapply(exp(r0*x),x)}));
 > simplify(%); # 0
 > rrcn:=normal(rrc);
 > simplify(denom(rrcn));
 > evalf(%); # 0.
 > simplify(numer(rrcn)); # 0

The problem appears to be Maple's ability to show that some quantity is actually zero.

Let us do a similar concrete test of the case a = 3  ( > 1/ln(2) ).

We didn't choose a = 2 because it is rather special: rm1 = -1 and r0 = 0 are found without evalf.

 > evalf(1/ln(2));
 > eval(rm1,a=3);
 > evalf(%);
 > rrc:=expand(eval[recurse]((lhs-rhs)(value(EQ0)),{a=3,Z=unapply(exp(rm1*x),x)}));
 > evalf(rrc); # Float(undefined)* ...
 > rrcn:=normal(rrc);
 > simplify(denom(rrcn));
 > evalf(%); # 0.
 > simplify(numer(rrcn)); #0

We look at r0.

 > eval(r0,a=3);
 > evalf(%);
 > rrc:=expand(eval[recurse]((lhs-rhs)(value(EQ0)),{a=3,Z=unapply(exp(r0*x),x)}));
 > simplify(%); #0
 > rrcn:=normal(rrc);
 > simplify(denom(rrcn));
 > evalf(%); # 9.295 ...
 > simplify(numer(rrcn));

Thus we (believe we) have found the following one parameter family of solutions to the inhomogeneous integral equation EQ:

 > SOL:=Zsol+C*piecewise(a<=0,0,a<=1/ln(2),exp(rm1*x),exp(r0*x));

A direct check:

 > rr:=expand(eval((lhs-rhs)(value(EQ)),Z=unapply(SOL,x))):
 > simplify(%); #  0

An animation in C for a = 1/2 on x = 1.5 .. 3.5:

 > SOL1:=eval(SOL,a=1/2);
 > evalf(SOL1);
 > plots:-animate(plot,[SOL1,x=1.5..3.5],C=-2000..2000,trace=24,caption=typeset("a = ",1/2));

An animation in C for a = 3 and on x = -1 .. 2:

 > SOL2:=eval(SOL,a=3);
 > evalf(SOL2);
 > plots:-animate(plot,[SOL2,x=-1..2],C=-1..1,trace=24,caption="a = 3");
 > ###################################################################################################################

We can single out one of these solutions by imposing an initial condition.

We shall assume that a > 0 from now on, since otherwise the solution SOL has no C:

 > eval(SOL,x=0) assuming a>0;

If Z(0) = 0 is required then C is given by C = C0, where:

 > C0:=solve(%,C);

The solution Z(x) = SOL is then

 > SOL0:=eval(SOL,C=C0) assuming a>0;

If furthermore a = 1/2 the solution is

 > SOL12:=eval(SOL0,a=1/2);
 > plot(SOL12,x=0..10);
 > limit(SOL12,x=infinity);
 > Digits:=15:
 > eps:=1e-5: plots:-animate(plot,[SOL0,x=0..10],a=eps..1/ln(2)-eps);

To find the value of SOL0 for a = 1 we must take the limit:

 > SOL1:=limit(SOL0,a=1);
 > limit(SOL1,x=infinity); # 3
 > plot([SOL1,3],x=0..20,linestyle=[1,3],color=[red,black],caption="a = 1");

The value of SOL0 at 1/ln(2) is found by taking the limit:

 > SOL_ln:=limit(SOL0,a=1/ln(2));
 > L_ln:=limit(SOL_ln,x=infinity); # 3*ln(2)/(-1 + 2*ln(2))
 > evalf(L_ln);
 > plot([SOL_ln,L_ln],x=0..20,linestyle=[1,3],color=[red,black],caption=typeset("a = ",1/ln(2)));
 >
 > ## No way...

I have corrected some syntactical errors. But dsolve cannot handle this analytically.

Certainly numerically, though:

```restart;
ode:=m*diff(x(t),t,t)= - m * A * sin(2*Pi*f*t) - piecewise( abs(x(t)) < x__max, 0, abs(x(t)) >= x__max, -k*(abs(x(t))-x__max)*signum(x(t))) - diff(x(t), t);
dsolve(ode); #NULL
nms:=indets(ode,And(name,Not(constant))) minus {x,t};
res:=dsolve({ode,x(0)=x0,D(x)(0)=x1},numeric,parameters=[nms[],x0,x1]);
res(parameters=[A=1,m=1,f=1,x__max=2,k=2,x0=0,x1=1]); # Example
plots:-odeplot(res,[t,x(t)],0..20);
``` ## Only rank 2...

Try this:

```restart;
Rok := Vector([2013, 2014, 2015, 2016, 2017, 2018], datatype = float);

TrzbyCelkemEmco := Vector([1028155, 1134120, 1004758, 929584, 995716, 1152042], datatype = float);
########
KubickaTrzby := Statistics:-PolynomialFit(3, Rok, TrzbyCelkemEmco, x);
p1:=plot(KubickaTrzby,x=2013..2018):
p2:=plot(Rok,TrzbyCelkemEmco,color=blue):
plots:-display(p1,p2); # Not impressive
########
pol:=a*x^3+b*x^2+c*x+d;
resid:=eval~(pol,x=~Rok)-TrzbyCelkemEmco; # Residuals
A,B:=LinearAlgebra:-GenerateMatrix(convert(resid,list),[a,b,c,d]);
op(1,A); # Dimension of A: 6x4
LinearAlgebra:-Rank(A); # 2 only
```

The rank could in principle have been min(6,4) = 4. That would have meant full rank.

## A start...

To get an exact answer to this question is certainly not easy to say the least.

But the beginnings of an attempt is this:

```restart;
fS := -D2*x1*x2-D1*x1-x1*x3+S;
fV := D2*x1*x2-D4*x2+x1*x3;
fC := -D6*x3*x4-D5*x3+x2;
fR := D6*x3*x4-x4;
eqs:=solve({fS,fV,fC,fR},{x1,x2,x3,x4}); # Notice RootOfs
L:=map(subs,[eqs],[x1,x2,x3,x4]); # The equilibrium points
### The Jacobian matrix as a function of x1,x2,x3,x4:
J:=unapply(VectorCalculus:-Jacobian([fS,fV,fC,fR],[x1,x2,x3,x4]),x1,x2,x3,x4);
J(x1,x2,x3,x4); # View
J(op(L)); # Point number 1
EV1:=LinearAlgebra:-Eigenvalues(%);
### The first two are negative since D1>0.
### But the two other depend on the actual size of your parameters as seen here:
### For these values the system is asymptotically stable:
eval(EV1,{D1=0.1,D2=1,D4=0.1,D5=.1,S=.0001});
eval(EV1,{D1=0.1,D2=1,D4=0.1,D5=.1,S=.0001});
### For these values the system is unstable:
eval(EV1,{D1=0.1,D2=1,D4=0.1,D5=.1,S=1});
eval(EV1,{D1=0.1,D2=1,D4=0.1,D5=.1,S=1});
```

This just illustrated the problem for point 1.

## dsolve...

Without specifying Q concretely you get:

```restart;
ode:=diff(v(x),x,x)+diff(v(x),x)*cot(x)-v(x)*(cot(x)^2+nu)=a^2*Q(x)/D1;
dsolve(ode);
```

Notice that there is an unevaluated integral involved simply because Q is not known.

## No solution...

That NULL is returned means that solve couldn't find a solution or that none exists.

In fact the latter is the case here with your parameters so far just names.
Your system is linear in the a-variables, the coefficient matrix A is 9x9 and its rank is only 8.
Thus for general right hand side B in A.x = B there is no solution for x = <a,..., a >.

```restart;

omega := v/h;
t := sum(a[j]*x^j, j = 0 .. 6)+a*cos(omega*x)+a*sin(omega*x);
r1 := diff(t, x\$2);
r2 := diff(t, x\$4);
c1 := eval(t, x = q+2*h) = y[n+2];
c2 := eval(r1, x = q) = f[n];
c3 := eval(r1, x = q+h) = f[n+1];
c4 := eval(r1, x = q+2*h) = f[n+2];
c5 := eval(r1, x = q+3*h) = f[n+3];
c6 := eval(r2, x = q) = g[n];
c7 := eval(r2, x = q+h) = g[n+1];
c8 := eval(r2, x = q+2*h) = g[n+2];
c9 := eval(r2, x = q+3*h) = g[n+3];
b1 := seq(a[i], i = 0 .. 8);
#####
A,B:=LinearAlgebra:-GenerateMatrix([c1, c2, c3, c4, c5, c6, c7, c8, c9],[b1]);
LinearAlgebra:-Dimension(A); #9x9
LinearAlgebra:-Rank(A); # 8
```

## Bug in solve: Use fsolve or workaround i...

Use fsolve instead. Your function is strictly increasing on the real axis, so the equation f(x) = 5 has at most one real solution.
It is odd and tends to infinity as x -> infinity. Thus there is exactly one real solution.

fsolve(f(x)=5,x);
############## Workaround ############

```restart;
f := x -> x*sqrt(4*x^2 + 1) + arcsinh(2*x);
convert(f(x),ln);
solve(%=5,x);
allvalues(%); # Only one
evalf(%);
#### With 5.0 instead of 5:
restart;
f := x -> x*sqrt(4*x^2 + 1) + arcsinh(2*x);
convert(f(x),ln);
solve(%=5. ,x);
```

## No solution...

You obviously by pi mean the mathematical constant, which in Maple is Pi.

Ideally, you would just do:

```restart;
ode := diff(y(x), x, x) + (n*Pi)^2*y(x) = A^3*sin(n*Pi*x)^3;

dsol1 := dsolve({ode, y(0) = 0, y(1) = 0}) assuming n::integer ;
```

It results in a division by zero error.

Are you sure that there is a twice differentiable solution?
Try letting n and A both be 1 for a concrete example:

```dsolve({eval(ode,{n=1,A=1}), y(0) = 0, y(1) = 0});
```

No solution is found.
Worse:

```ode1:=eval(ode,{n=1,A=1});
sol:=dsolve({ode1, y(0) = 0, D(y)(0) = y1});
eq:=eval(rhs(sol),x=1);
```

Since eq = 3/(8*Pi) is independent of y1, no value of y'(0) = y1 will make y(1) = 0.

Thus there is no solution to your problem (here only shown for n=1 and A=1).

A general right hand side for n=1:

```restart;
ode := diff(y(x), x, x) + Pi^2*y(x) = f(x);
sol:=dsolve({ode,y(0)=0,D(y)(0)=y1});
eq1:=eval(sol,x=1);
## Examples:
value(eval(eq1,f=sin)); # not zero
eval(eq1,f=sin^2);
value(%);
evalf(%);
```

You need f to satisfy: int(sin(Pi*x)*f(x), x = 0 .. 1) = 0.

A simple example of a left hand side for which there is a solution is f(x) = x - 1/2, because eq1 is satisfied:

```ode1:=eval(ode,f(x) = x-1/2);
dsolve({ode1,y(0)=0,y(1)=0});
```

You will notice that now there are infinitely many solutions.

## solve, fsolve...

You have a polynomial system with 7 polynomials in the variables indets(sys,name) = {omega, x, a[-1], a, a, b[-1], b}.

The coefficients are just integers. So in principle you can just try
solve(sys);
that will be understood by Maple as solve(sys,{omega, x, a[-1], a, a, b[-1], b} );
Does the computation ever come to an end?
Try fsolve the same way, but be aware that you are only getting one solution out of many
fsolve(sys);
fsolve(sys, complex);

## You could also try SolveTools:-PolynomialSystem:

SolveTools:-PolynomialSystem(sys,{omega, x, a[-1], a, a, b[-1], b},maxsols=20);

## Syntactical problems cleared...

The following may not do what you intended, but at least the syntactical problems are gone.

You don't need any assumption on a.

```restart;
N:=10;
b:=1: b:=-1: b:=3: #Example
####
for n from 3 to N do
if (n mod 3) = 0  then
b[n] :=b
elif   (n mod 3) = 1 then
b[n]:=b
elif  (n mod 3) = 2  then
b[n]:=b
end if
end do;
####
seq(b[i],i=1..10);
####
u:=1: #Example
####
## You may want to assign a value to a before this, but you don't have to:
for n from 1 to N do  u[n+1]:=expand(a*u[n])+b[n]  end do;
### Seeing what you have for a = 0.4:
seq(eval(u,a=0.4)[n],n=1..N);
```

## Elementwise or map or directly...

Beginning with Maple 2018, you can just do diff(A, x).
In earlier versions, but later than Maple 12 you can use diff elementwise: diff~(A,x).
In Maple 12 and earlier you can use map: map(diff, A, x).

map works in all versions. Elementwise works in all versions from Maple 13 and up.
Example:

```A:=Matrix( [[sin(x),4*exp(3*x)],[cos(5*x),sinh(cos(x))]] );
diff(A,x);        # Maple 2018 and later
diff~(A,x);      # Maple 13 and later
map(diff,A,x); # All versions
```

The error you get in (e.g.) Maple 2017 from trying diff(A,x) is

Error, non-algebraic expressions cannot be differentiated

The error you get in Maple 12 from doing diff~(A,x) is more obscure:

Error, missing operator or `;`

## value...

Use value. It turns the inert Int into int.

Example:

```V:=Vector[row]([seq(Int(x^n/L^(n+2),x=0..L),n=1..5)]);
value(V);
```

You are using evalm. That was used in the old and deprecated linalg package.

To get a matrix similar to yours you can do

```V^+.V;
value(%);```

## Using Maple 12...

Since you are using Maple 13, I tried in Maple 12 (don't have 13):
My guess is that there is not much difference between the two.

```## Maple 12.
restart;
eqs:={(13/4)*m-(7/4)*n-3 = 0,-(17/2)*n*2^n +34*m= 0};
plots:-implicitplot(eqs,m=-1..3,n=-4..3);
sol:=solve(eqs);
evalf(sol); # Your desired solution, but as floats.
A:=allvalues(sol);
evalf(A); # Misses the other!
E:=eliminate(eqs,m);
eq_n:=op(E);
plot(eq_n,n=-3..3); # Looks like 2 solutions (at least)
diff(eq_n,n);
z_n:=solve(%,n);
#The derivative has one zero only.
#Thus there cannot be more than 2 (real) solutions for n (and hence for (m,n) ).
evalf(z_n);
plot(diff(eq_n,n),n=-1..1);
fsolve(eq_n);
fsolve(eq_n,n=1);
```

So there are two real solutions for (m,n).

## List of lists...

The initial conditions should be given as a list of lists even when there is only one point ( here (1,1) ).

```restart;
with(DEtools):
LC := [diff(x(t), t) = x(t)*(1-x(t)^2-y(t)^2)+y(t)*((1-x(t))^2+y(t)^2), diff(y(t), t) = y(t)*(1-x(t)^2-y(t)^2)-y(t)*((1-x(t))^2+y(t)^2)];
p1 := DEplot(LC, {x(t), y(t)}, t = 0 .. 100, [[x(0) = 1, y(0) = 1]]); ## list of list
```

﻿