Preben Alsholm

13728 Reputation

22 Badges

20 years, 248 days

MaplePrimes Activity


These are answers submitted by Preben Alsholm

There are several equality signs that (I suppose) ought to be assignments (:=) instead.
In the assignment to PositiveEQU2, RT should be R*T.
Besides that, using R and R[f] (and more like that) is very likely to cause serious problems.
To see that you may try:
R := 8.3145; F := 96485.3365; T := 298.15; R[f] := 0;
whattype(eval(R));
eval(R);
I suggest changing the first name to Rg or similar. Using R[f] and the other indexed variables ought to be safe.
However, even after these changes, try
indets(BatSys,name);
I found the output to contain besides t and x also R[pN], R[pP].
R[pN] and R[pP] are surprising as they seem to be assigned numerical values.
However, try:
lprint(NegativePDE2);
and compare with
lprint(NegativePDE1);
The problem is with the 2d input. I never use it.
Certainly indets(BatSys,name) should only return {t, x}. So before that is the case, pdsolve/numeric won't work.

The procedure works if you change it to read:
p:=proc(pp2) local res,F0,F1,F2;
if not type([pp2],list(numeric)) then return 'procname(_passed)' end if:
res := dsolve({eq1=0,subs(p2=pp2,eq2)=0,eq3=0,u(0)=0,u(1-zet)=0,phi(0)=0.008,T(0)=0,D(T)(0)=1}, numeric,output=listprocedure);
F0,F1,F2:=op(subs(res,[u(eta),phi(eta),T(eta)]));
evalf(2/(1-zet^2)*Int((1-eta)*(F1(eta)*rhop+(1-F1(eta))*rhobf)*( F1(eta)*rhop*cp+(1-F1(eta))*rhobf*cbf )/(F1(eta)*rhop+(1-F1(eta))*rhobf)*F0(eta),eta=0..1-zet))-pp2
end proc;

I added output=listprocedure, and removed the line res(parameters=[pp2]) as that doesn't belong here when the parameters option is not (and cannot be) used.

m:=4:
A:=LinearAlgebra:-BandMatrix([0$m,1,seq(a[i],i=1..m)], m,m+1 );
B:=Vector(m+1,i->b[i-1]);

Since K is allowed to depend on the a's, there are quite a lot of solutions:
restart;
with(LinearAlgebra):
K:=Matrix(3,symbol=b);
V:=Vector(3,symbol=a);
eq1:=a[1]+a[2]*a[3]^2*a[1]+a[2]*a[3]-a[1]^2;
eq2:=a[3]+a[1]^3+a[2]*a[3]+a[2]^2-a[1]^2;
eq3:=a[1]^3+a[2]^3+a[1]*a[2]*a[3]+a[1]*a[2]-a[3];
W:=<eq1,eq2,eq3>;
K.V=W; #Need to solve for K, i.e. the b's:
Equate(K.V,W);
res:=solve(%,{seq(seq(b[i,j],j=1..3),i=1..3)}); Using solve
K1:=subs(res,K);
K1.V;

The first example I don't understand: Where is the dependence on t? Are x and y supposed to be dependent on t?

The second equation is a partial differential equation and may also be written:

pde:=Diff(f(x,y), x,x,y,y)*x+Diff(f(x,y), x,y)*y = 0;
#The solver here is pdsolve:
pdsolve(pde);

Doesn't it follow from the boundary conditions at x = 0 and from the PDE (without solving it) that diff(u(x,t),x,x) and diff(v(x,t),x,x) are both zero at x = 0 for all t?
However,
Robert Lopez gave an idea in the link given at the top. Introduce two new functions u2x and v2x:

PDE2 := {diff(u(x,t),t)=-1/20*v2x(x,t),diff(v(x,t),x,x)=v2x(x,t),
diff(v(x,t),t)=u2x(x,t),diff(u(x,t),x,x)=u2x(x,t)};
pds2 := pdsolve(PDE2, IBC, numeric, spacestep=1/50);
VAL:=pds2:-value(output=listprocedure);
VAL(0.5,3);
u2xVal,v2xVal:=op(subs(VAL,[u2x(x,t),v2x(x,t)]));
u2xVal(0,5);
v2xVal(0,1.2345);
plot3d(u2xVal,0..1,0..5);
plot3d(v2xVal,0..1,0..5);

In view of the fact that roots of polynomials are often nasty looking when they can be found symbolically at all, it may be worth your while rewriting the system as a first order linear system X' = A.X and use linear algebra. The solution is X = exp(A*t).X0, where X0 is  the initial vector X(0) = X0.
The point here is that if you replace A with its floating point version (i.e. use evalf(A) instead of A) then you are sure that eigenvectors and eigenvalues can be found so that MatrixExponential wil be able to return a result.
A similar device will not work when using dsolve since it turns floating point coefficients into rationals before doing any serious work.

restart;
eq1 := diff(x1(t),t,t) = -2*x1(t)+x2(t)-diff(x1(t),t) ;
eq2 := diff(x2(t),t,t) = -x2(t)+x1(t);
ics:= x1(0) = 1, x2(0) = 0, D(x1)(0) = 0, D(x2)(0) = 0;
#Turning the system into a first order system:
DEtools[convertsys]([eq1,eq2],[ics],[x1,x2],t,Y,YP);
subs(%[1],[YP[i]$i=1..4]);
A:=LinearAlgebra:-GenerateMatrix(%,[Y[i]$i=1..4])[1];
X:=<x1(t),x1p(t),x2(t),x2p(t)>; #The naming follows the output from convertsys
diff~(X,t) = A.X; #The system (just for show)
LinearAlgebra:-MatrixExponential(evalf(A),t).<1,0,0,0>; #Notice evalf(A) instead of A
res:=fnormal~(%); #Removing roundoff errors
plot([res[1],res[3]],t=0..30); #Plotting x1 and x2

I'm assuming that you have a couple of typos: k should have been k1 and in eq2 diff(x1(t),t,t) should have been diff(x2(t),t,t).
restart;
#k replaced by k1
eq1 := m*diff(x1(t), t, t)-f = 0; f := -k1*x1(t)+k2*(x2(t)-x1(t))-lambda*diff(x1(t), t);
#diff(x2(t), t, t) instead of diff(x1(t), t, t):
eq2 := m*diff(x2(t), t, t)-g = 0; g := -k2*(x2(t)-x1(t));
fcns := x1(t), x2(t);
ICS := {x1(0) = 1, x2(0) = 0, (D(x1))(0) = 0, (D(x2))(0) = 0};
sys := {eq1, eq2};
m := 1; k2 := 1; k1 := 1; lambda := 1;
sysdiff := sys union ICS;
#dsolve(sysdiff); #No luck: Arbitrary constants appearing
Sol_N := dsolve(sysdiff,{fcns},method=laplace);
res:=allvalues(value(Sol_N));
plot(subs(res,[x1(t),x2(t)]),t=0..30);



Maybe this small example using the parameters option can help.


restart;
eq:=diff(f(t),t,t)+f(t)=t/(t+p);
sol:=dsolve({eq,f(0)=1,D(f)(0)=0},numeric,parameters=[p],events=[[f(t)-1.2,halt],[f(t)+0.2,halt]]):
P:=proc(p) local res,q;
    sol(parameters=[p]);
    res:=op(subs(sol(15),[t,f(t),diff(f(t),t)]));
    q:=sol(eventstop);
    [p,res,q];
end proc:
P(1.4);
P(2.7);
P(10);

If the expression is not defined as a function (as in Kitonum's answer), then you can use diff and eval:

u:=x^2+2*x-3;
eval(diff(u,x),x=2);

You can still use diff and eval in the function case, but it is not as elegant as using D.
Continuing from u:

f:=unapply(u,x);
eval(diff(f(x),x),x=2);
#Compare with the simpler:
D(f)(2);

The problem is that you have a multivariate polynomial with floating point coefficients.
See the help page for factor:
?factor
For a univariate polynomial there is no problem:
u1:=1.2*x+3.4*x^2;
factor(u1);
#Now a polynomial in the variables x and a:
u:=1.2*a*x+3.4*x^2;
factor(u); #Doesn't factor
#Turning the floats into rationals:
factor(convert(u,rational,exact)); #Factors
evalf(%);
#Turning the floats into symbols instead:
evalindets(u,float,x->convert(x,symbol));
factor(%);
evalindets(%,symbol,parse);

The error message tells the story. So you are on your own.
Here is a very simple example which may or may not help you in your more complicated situation.

ode:=diff(x(t),t,t,t)=sin(x(t)); #Just one equation
dsolve({ode,x(0)=1,x(1)=0.5,x(2)=2},numeric); #Error as expected
#Now only 2 points, but with an additional requirement at 0, in this case D(x)(0) = a, where a is to be determined so that x(1) = 0.5:
p:=a->dsolve({ode,x(0)=1,x(2)=2,D(x)(0)=a},numeric,output=listprocedure);
p1:=proc(a) local X1; X1:=subs(p(a),x(t)); X1(1)-0.5 end proc;
#Just exploring the ground:
plots:-odeplot(p(-1),[t,x(t)],0..2);
#Now looking for a zero for p1, i.e. a value for a for which x(1) = 0.5 using a = -1 as an initial guess:
A:=fsolve(p1,-1);
p(A)(1);
plots:-odeplot(p(A),[t,x(t)],0..2);

If you move one line and modify it, I think you got it:
for A  from .1 by .01 to 3 do  
   for B  from .1 by .01 to 3 do
      l1:=(exp(B)-A)-1;
      l2:=(exp(B)-A/2)-1.5;
      if l1>0 and l2>0 then print(A,B); break end if
   end do
end do;

#Or if you don't want to copy from the screen:

S:=NULL:
for A  from .1 by .01 to 3 do  
   for B  from .1 by .01 to 3 do
      l1:=(exp(B)-A)-1;
      l2:=(exp(B)-A/2)-1.5;
      if l1>0 and l2>0 then S:=S,[A,B]; break end if
   end do
end do;
S; #A sequence of pairs [A,B]

As far as dsolve is concerned your syntax is not quite right. Try this:

restart;
mystery := proc (z) evalf((1+I)*z) end proc;
F1 := proc (N, t, Y, YP) local tmp;
   tmp:=mystery(Y[1]+I*Y[2]);
   YP[1] := Re(tmp);
   YP[2] := -Im(tmp)
end proc;
curve := dsolve(numeric,procedure = F1, number = 2, procvars = [x(t), y(t)], start = 0, initial = Array([-.2,.3]));
plots:-odeplot(curve,[x(t),y(t)],-1..1);
curve(1.2345);
#The last output with the .0*I terms hints at the possible problem DEplot is having.
#At the moment I only have access to Maple 12, but there it seems that the parameters option doesn't work when the parameters don't appear in the system, in this case in F1.

dsys := {diff(x(t),t)=y(t),diff(y(t),t)=-x(t),x(0)=1,y(0)=0};

dsn := dsolve(dsys,numeric,output=listprocedure);

X,Y:=op(subs(dsn,[x(t),y(t)]));
evalf(Int(X(t),t=0..1));   
#For your simple system you could have found the integral from the second equation much simpler:
Y(0)-Y(1);

If the system is more complicated you can still use the latter method simply by introducing z(t) by
diff(z(t),t) = x(t) with z(0) = 0. Then z(1) = int(x(t),t=0..1):

restart;
dsys2 := {diff(x(t),t)=y(t),diff(y(t),t)=-sin(x(t)),diff(z(t),t)=x(t),x(0)=1,y(0)=0,z(0)=0};
dsn2 := dsolve(dsys2,numeric,output=listprocedure);
X,Y,Z:=op(subs(dsn2,[x(t),y(t),z(t)]));
Z(1);
#Compare with:
dsys := {diff(x(t),t)=y(t),diff(y(t),t)=-sin(x(t)),x(0)=1,y(0)=0};
dsn := dsolve(dsys,numeric,output=listprocedure);
X,Y:=op(subs(dsn,[x(t),y(t)]));
evalf(Int(X(t),t=0..1));   





First 98 99 100 101 102 103 104 Last Page 100 of 160