Preben Alsholm

12103 Reputation

22 Badges

16 years, 195 days

MaplePrimes Activity


These are answers submitted by Preben Alsholm

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

 

At least in your simple example converting to Heaviside works:
 

F := a*X+piecewise(Y<0, X, b);
FH:=convert(F,Heaviside);
type(FH, linear(X));

 

There is no need for that kind of requirement.

Your ode has two constant solutions: x(t) = a and x(t) = b/2.

Your ode satisfies requirements for existence and uniqueness of intitial value problems.
Thus a solution starting with x(0) = 0 cannot ever hit x(t) = a or b/2, since this would violate uniqueness.

I once asked Edgardo Cheb-Terrab (in this forum, I think) why the unknown function (here y(t)) had to be given when the series and laplace methods are used while it is not in general necessary otherwise.
In his answer he explained that the code for those two methods were not written by him, but by other people.

So probably the simple answer to your question is that those "other people" didn't implement the eval-variant.
If you try:

dsolve([ode,bc_form_1],y(t),method=laplace);

then you get a rather misleading error message:

Error, (in dsolve/inttrans/solveit) transform method can only be applied to linear ODE
 

 

 

 

Actually Maple 2018 comes up with the correct solution , {T(x, t) = x}, {V(x, t) = 1+cos(w*t)} ,

but Maple 2019 says no solution (i.e. NULL):

restart;
sys_Pde := diff(V(x, t), x, x) = 0, diff(T(x, t), x, x) = 0; 
BC := eval(diff(V(x, t), x), x = 1) = 0, eval(V(x, t), x = 0) = 1+cos(w*t), eval(T(x, t), x = 0) = 0, eval(T(x, t), x = 1) = 1;

res:=pdsolve({sys_Pde}); # General solution
pdsolve({BC, sys_Pde}); # 2018 OK, 2019 NOT OK.
## Check:
Vsol:=subs(res,V(x,t));
Tsol:=subs(res,T(x,t));
res1:=diff(Vsol,x)=0; # So _F3(t) = 0.
res2:=eval(Vsol,x=0)=1+cos(w*t); # So _F4(t) = 1 + cos(w*t).
res3:=eval(Tsol,x=0)=0; # So _F2(t)=0.
res4:=eval(Tsol,x=1)=1; # So _F1(t)+_F2(t)=0. Thus _F1(t) = 1.  (!)
sol:=eval(res,{res1,res2,res3,res4});  # Both missing that _F1(t)=1.
pdetest(sol,{BC, sys_Pde});                # Both missing that _F1(t)=1       

Which Maple version are you using?

OK, I see that the error that you report comes up in Maple 2017.

Notice that your procedure gra only works for the integers 1..200000 (nmax). It doesn't work for e.g. the float 88.0.
The easiest way out of that is to use plots/insequence:

plots:-display([seq(gra(h),h=0..200000,1000)],insequence=true);

 

I shall assume that E = R^3, so just do:
 

int(exp(-x^2-y^2/4-z^2/9),[x=-infinity..infinity,y=-infinity..infinity,z=-infinity..infinity]);

Answer: 6*Pi^(3/2)

If E is not R^3 you should tell us what it is.

5 6 7 8 9 10 11 Last Page 7 of 151