Preben Alsholm

11749 Reputation

22 Badges

15 years, 295 days

MaplePrimes Activity


These are answers submitted by Preben Alsholm

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);

 

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:

 

 

 

https://www.mapleprimes.com/questions/228874-Collocation-Method-For-Integral-Equation-#answer266560

Preben Alsholm, March 17 2020.
Additions made at the bottom March 18 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[15](%);

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[15](%);

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[15](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[15](SOL2);

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

###################################################################################################################

Added March 18 2020:

 

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)));

 

 

 

 

 

Download IntegralEquations_delay_Z.mw

 

 

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);

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.

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;
### Your right hand sides:
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[1])); # 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[4],{D1=0.1,D2=1,D4=0.1,D5=.1,S=.0001});
eval(EV1[3],{D1=0.1,D2=1,D4=0.1,D5=.1,S=.0001});
### For these values the system is unstable:
eval(EV1[4],{D1=0.1,D2=1,D4=0.1,D5=.1,S=1});
eval(EV1[3],{D1=0.1,D2=1,D4=0.1,D5=.1,S=1});

This just illustrated the problem for point 1.

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.

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[0],..., a[8] >.
 

restart;

omega := v/h;
t := sum(a[j]*x^j, j = 0 .. 6)+a[7]*cos(omega*x)+a[8]*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

 

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);

 

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.

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

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[0], a[1], b[-1], b[0]} );
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[0], a[1], b[-1], b[0]},maxsols=20);

 

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[0]:=1: b[1]:=-1: b[2]:=3: #Example
####
for n from 3 to N do    
  if (n mod 3) = 0  then
    b[n] :=b[0]   
  elif   (n mod 3) = 1 then
    b[n]:=b[1]   
  elif  (n mod 3) = 2  then
    b[n]:=b[2]   
  end if  
end do;
####
seq(b[i],i=1..10);
####
u[1]:=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);

 

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 `;`

 

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(%);

 

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[2]);
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).

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

 

2 3 4 5 6 7 8 Last Page 4 of 148