Preben Alsholm

13471 Reputation

22 Badges

20 years, 298 days

MaplePrimes Activity


These are answers submitted by Preben Alsholm

Both of these work:

readdata("F:/test.txt",[string,integer,integer,integer]);
evalindets(%,string,()->NULL);
#ImportMatrix may be easier to use:
ImportMatrix("F:/test.txt",delimiter=" ");
convert(%,listlist);
evalindets(%,name,()->NULL);


Since the "correct" answer is supposed to be 3 lists M1, M2, and M3, you really have 9 real variables.

It seems to me that you have 3 equations that are sums of dotproducts, the first one being

rho2.M1 - rho1.M2 = 0.

Now working with vectors and matrices instead of lists you could do as follows (if I understand your intention):

restart;
F := [x1/(1+t^2)+x2, -(x1+x2)*(1+1/(1+t^2)), -x3];
#Let the rho's be row vectors:
rho1 := <diff(F[1],x1)|diff(F[2],x1)|diff(F[3],x1)>;
rho2 := <diff(F[1],x2)|diff(F[2],x2)|diff(F[3],x2)>;
rho3 := <diff(F[1],x3)|diff(F[2],x3)|diff(F[3],x3)>;
theta := Matrix([[rho2, -rho1, <0|0|0>],[rho3, <0|0|0>, -rho1],[<0|0|0>, rho3,-rho2]]);
#The proposed solution is checked:
Mx:=<-1,-1,0,1,0,0,0,0,-1>; #The 3 solutions stacked on top of each other.
theta.Mx;
#Well, indeed a solution, but there are infinitely many more:
LinearAlgebra:-LinearSolve(theta,<0,0,0>,free=f);

I tried with success

restart;
evalf(Sum(Zeta(p-2.)*2.9^p/GAMMA(.5*p+1.), p = 1000..infinity));
evalf(Sum(Zeta(p-2.)*2.9^p/GAMMA(.5*p+1.), p = 4 .. 1000));
restart;
evalf(Sum(Zeta(p-2.)*8^p/GAMMA(.5*p+1.), p = 1000..infinity));
evalf(Sum(Zeta(p-2.)*8^p/GAMMA(.5*p+1.), p = 4 .. 1000));

But don't know why Maple doesn't succeed with the entire range at once.
#Comment: Actually the following works:
restart;
evalf[250](Sum(Zeta(p-2.)*2.9^p/GAMMA(.5*p+1.), p = 4 .. infinity));
evalf(%);
#But this doesn't. It returns Float(infinity) (and takes a long time)
evalf[2000](Sum(Zeta(p-2)*8^p/GAMMA(.5*p+1.), p = 4 .. infinity));


convert(par,listlist);
map( ((x, y)-> x^2+y^2)@op , %);

What you have in the first integral is essentially this:

Int(x^j*exp(-a*x),x=0..infinity);
value(%) assuming a>0,j>-1;
simplify(%) assuming j::posint; #If that is indeed the case?
convert(%,factorial);



Try this:

eq1,eq2:=diff(f1(x),x)=0.25*b*(a/b-f1(x))^2, diff(f2(x),x)=0.25*b*(a/b-f2(x)+f1(x))^2;
ic1,ic2:=f1(0)=0,f2(0)=0;
dsolve({eq1,eq2,ic1,ic2});
simplify(%);
#The (escaped local) names _Z1~ and   _Z3~ (or similar) represent arbitrary integers:
indets(%,`local`);
If you put them equal to zero you get the same result as when you solve in two steps as done afterwards:
res1:=eval(%%,%=~0);
dsolve({eq1,ic1});
subs(%,{eq2,ic2});
res2:=dsolve(%);
subs(res1,f2(x))-subs(res2,f2(x));
simplify(%);
#In Maple 16 it simplifies to 0.


#Comment: If you do instead res1:=eval(%%,%=~1); you get a result for which f2(0) <> 0. Since the initial value has a unique solution we know in advance that either the solution doesn't in fact depend on the escaped locals (the _Z's) or the solution is wrong for some sets of _Z's.

You need to replace G by G(t,z[1],z[2],z[3]), but I don't see much of a chance for obtaining a solution with Maple:

Eq1 := subs(  G=G(t,z[1],z[2],z[3]),'((z[1]-1)*lambda+(z[3]-1)*gamma)*G+(1-z[1])*delta*(diff(G, z[1]))+(N*z[3]-z[2])*kappa*(diff(G, z[2]))+(1-z[3])*mu*(diff(G, z[3]))+(z[2]-z[1]*z[3])*beta*(diff(diff(G, z[3]), z[1])) = diff(G, t)' );
#pdsolve separates variables, but that is it:
pdsolve(Eq1);

Try this:

A:=Int(Int((1-x)*(1-y)*(x+y)*(x^n_1)*(y^n_2)*(1-x-y)^n_3, y=0..1-x),x=0..1);
#I have used Int instead of int so we see the integral before proceeding with a calculation.
#Decent looking results (here in Maple 12) are only available for concrete values of the exponents:
eval(A,{n_1=1,n_2=2,n_3=3});
value(%);
#But you may try
value(A) assuming posint;

Not many corrections were needed.
restart;
#A, B = Payoff matrices for players
A := Matrix(2, 2, [4, 2, 1, 3]); B := Matrix(2, 2, [1, 4, 3, 2]);
#T defines the time interval [0,T], s is the step length
T := 80; s := 1/5; k := [seq(i*s, i = 1 .. T)];
#Trajectories for different initial conditions
S1 := dsolve({z1(0) = .5, z2(0) = .18, diff(z1(t), t) = z1(t)*(1-z1(t))*(A[1, 2]-A[2, 2]-(A[1, 2]+A[2, 1]-A[1, 1]-A[2, 2])*z2(t)), diff(z2(t), t) = z2(t)*(1-z2(t))*(B[1, 2]-B[2, 2]-(B[1, 2]+B[2, 1]-B[1, 1]-B[2, 2])*z1(t))}, numeric,output=Array(k));
S2 := dsolve({z1(0) = .5, z2(0) = .1, diff(z1(t), t) = z1(t)*(1-z1(t))*(A[1, 2]-A[2, 2]-(A[1, 2]+A[2, 1]-A[1, 1]-A[2, 2])*z2(t)), diff(z2(t), t) = z2(t)*(1-z2(t))*(B[1, 2]-B[2, 2]-(B[1, 2]+B[2, 1]-B[1, 1]-B[2, 2])*z1(t))}, numeric,output=Array(k));
#Trajectories
g1 := [seq([S1[2, 1][i, 3], S1[2, 1][i, 2]], i = 1 .. T)]:
g2 := [seq([S2[2, 1][i, 3], S2[2, 1][i, 2]], i = 1 .. T)]:
plot([g1, g2], axes = boxed, color = [black,blue]);

I don't see any immediate solution, but the following at least works.


C:=Matrix([[D[1],0,D[2]],[0,D[2],D[1]]]);
B:=;
freeze~(B);
C.%;
evalindets(%,`*`,apply@op);
thaw~(%);
%;

#Putting this into a procedure:

`&p`:=proc(d,v) local a,b;
   a:=d.freeze~(v);
   b:=thaw~(evalindets(a,`*`,apply@op));
   eval(b)
end proc;
   
C &p B;

#####Later addition:
In the more general case where your differential operators in the matrix are polynomials in D[1] and D[2] you will have to freeze C too:

Take as a simple example the matrix A below which has the linear differential operator T as one of its elements:

T:=-3*D[1]+7*D[2]+(f->6*f);
#Testing it on a function g:
T(g);
A:=Matrix([[T,0,D[2]],[0,D[2],D[1]]]);
#Now just modify the procedure above so that A will be frozen too:
`&p`:=proc(d,v) local a,b;
   a:=freeze~(d).freeze~(v);
   b:=thaw~(evalindets(a,`*`,apply@op));
   eval(b)
end proc;
   
A &p B;


I have corrected typos and whatever to get the following version, which works.

restart:

Eq1:=diff(f(eta),eta,eta,eta)+0.5*f(eta)*diff(f(eta),eta,eta)=0;

Eq2:=diff(theta(eta),eta,eta)+0.5*Pr*f(eta)*diff(theta(eta),eta)=0;

bc1:=f(0)=0, D(f)(0)=0, D(f)(10)=1;

#No []
fixedparameter:=Pr=0.72;

#eval not evalf
Eq3:=eval(Eq2,fixedparameter);

#a*... ? and theta(eta) probably should be theta(0) or ?
bc2:=theta(10)=0, D(theta)(0)= -a*(1-theta(0));

L:=[0.05, 0.1, 0.2, 0.4, 0.6, 0.8, 1, 5, 10, 20];

for k from 1 to 10 do

#bc1,bc2 (no s)
R:=dsolve(eval({Eq1, Eq3, bc1, bc2},

a= L[k]), [f(eta), theta(eta)], numeric, output = listprocedure);


Y||k :=rhs(R[5]); Yp||k := -rhs(R[6]);

end do:

#plot, not plots
plot([Y||(1..10)],0..6,labels= [eta, thet(eta)]);

[Yp||(1...10)(0)];
[Y||(1..10)(0)];

(Comment: I changed the names of the procedure for my own use. That of course is irrelevant.)
Using the result from LSSolve with Le = 8 as an initialpoint for LSSolve when Le=10 I was able to get a result:

With Le=5 first do

res:=Optimization:-LSSolve([pPhi,pTheta,pF1,pS]);
inipt:=convert(res[2],list);
#Keep inipt for later use
#Next set Le = 8 and do
Optimization:-LSSolve([pPhi,pTheta,pF1,pS],initialpoint=inipt);
#Again keep the result for later use in inipt2, say.
#Finally set Le = 10 and do
Optimization:-LSSolve([pPhi,pTheta,pF1,pS],initialpoint=inipt2);

#You can get a result for Le = 8 without option 'initialpoint', but that requires adding the option maxfun=0 in the call to dsolve/numeric. It takes a while then.

 

#Added July 5:
Here is a faster version which uses option remember in a procedure 'paramset', whose only job is to set the parameters.
Since only the procedural version of LSSolve is used there is no need for the lines about returning unevaluated on nonnumeric input.

restart;
eq1:=diff(f(x), x$3)+diff(f(x), x$2)*(A*f(x)+x/2)+diff(f(x), x$1)*(1- A*diff(f(x),x$1))+A-1  ;
eq2:=diff(s(x), x$2)+diff(s(x), x$1)*(A*f(x)+x/2)+s(x)*(1-2*diff(f(x),x$1));
eq3:=diff(theta(x), x$2)/Pr+A*( f(x)*diff(theta(x),x$1))+x*0.5*diff(theta(x), x$1)+Nb*diff(theta(x),x$1)*diff(phi(x),x$1)+Nt*diff(theta(x),x$1)^2;
eq4:=diff(phi(x), x$2)+(A*f(x)+x*0.5)*diff(phi(x), x$1)*Le+Nt*diff(theta(x),x$2)/Nb;

b := 8:  A:=2: Nb:=0.1: Nt:=0.1: Pr:=1:
#I let 'Le' be unassigned, but include it as a parameter in dsolve.
#Le:=8:
solpar:=dsolve({eq1=0,eq2=0,eq3=0,eq4=0, f(0)=0, D(f)(0)=0, D(D(f))(0)=f2,s(0)=1, D(s)(0)=s1,theta(0)=1, D(theta)(0)=t1,phi(0)=1, D(phi)(0)=p1}, parameters=[Le,f2,t1,p1,s1],numeric,output=listprocedure,maxfun=0):
phin,thetan,sn,f1n:=op(subs(solpar,[phi(x),theta(x),s(x),diff(f(x),x)])):

paramset:=proc(Le,f2,t1,c1,s1) option remember; solpar(parameters=[Le,f2,t1,c1,s1]); NULL end proc;
#pPhi computes phi(b)
pPhi:=proc(f2,t1,p1,s1)
paramset(Le,_passed);
phin(b)
end proc:
#pTheta computes theta(b)
pTheta:=proc(f2,t1,p1,s1)
paramset(Le,_passed);
thetan(b)
end proc:

#pS computes sn(b)
pS:=proc(f2,t1,p1,s1)
paramset(Le,_passed);
sn(b)
end proc:

#pF1 computes f1n(b)-1
pF1:=proc(f2,t1,p1,s1)
paramset(Le,_passed);
f1n(b)-1
end proc:

forget(paramset); #Erases the memory table just in case you do this several times
t0:=time():
inipt:=table():
inipt[4]:=[]:
for Le from 5 to 10 do
  res:=Optimization:-LSSolve([pPhi,pTheta,pF1,pS],initialpoint=inipt[Le-1]);
  inipt[Le]:=convert(res[2],list)
end do;
time()-t0;

#There is an interesting problem when Le=9.

Try

indets({Solseq1},`local`);
indets({Solseq2},`local`);

The copied version will after execution refer to the global names d and `a[1]`.

restart;
f:=(-31250000 + 312500 *(x + y) +
 625 *(4* a* x - 4* x^2 - 5* x* y + 4 *(b - y)* y) +
 25* x* (-a* b + b* (x - y) + y^2))/((-500 + 5* a)* (500 - 5* b));
g:=(-1/(1200* (100 - a)* (-100 + b)))*(-2 *a^3* (-100 + b) + 200* b^3 - b^4 +
   6* a* (-100 + b)* (200 - x)* x + 120000* x^2 - 800* x^3 + x^4 +
   120000* y^2 - 1200* x* y^2 - 800* y^3 + 8 *x* y^3 +
   6* b* (x^3 - 20000* y - x* y^2 + 100* (-2* x^2 + 2* x* y + y^2)));
#I hope I got your f and g right.
pxy:=proc(aa,bb) global f,x,y; local F,res;
F:=subs(a=aa,b=bb,f);
res:=Optimization:-Minimize(F,x=aa..100,y=bb..100);
op(subs(res[2],[x,y]))
end proc;
#Testing:
pxy(.4,5);
G:=unapply(g,a,b,x,y);
gg:=(a,b)-> G(a,b,pxy(a,b));
Testing:
gg(1,2);
Optimization:-Maximize(gg,0..100,0..100);

Does the following simple example suit your situation. The original equation is simply diff(x(t),t)=(-x(t)^2-1).

sys:=diff(x(t),t)=(-x(t)^2-1)*piecewise(x(t)>=0,1,0);
sol:=dsolve({sys,x(0)=1},numeric);
plots:-odeplot(sol,[t,x(t)],0..10,thickness=2);

First 120 121 122 123 124 125 126 Last Page 122 of 158