Preben Alsholm

MaplePrimes Activity


These are answers submitted by Preben Alsholm

##TRY:
combinat:-choose([$(0..n-1)], k);
add(mul~(combinat:-choose([$(0..n-1)], k))) ;

You aren't going to get explicit solutions for this equation. But you do get some good information from an initial symbolic treatment.
After that use numerical solution.
Notice that because a'(t) is squared dsolve will first have to find the 2 solutions for a'(t) and then solve those, in this case using separation of variables. The first result that pops up besides the 2 mentioned corresponds to the exceptional case when separating variables: that a denominator can be zero.
 

restart;
ode:=3*a(t)-a(t)^2-3*z*diff(a(t), t)^2+c*a(t)+w/a(t) = 0; # c*k replaced by just c
dsolve(ode); # 3 results
value([%]); #Real large so forget it.
sol:=dsolve(ode,implicit); # Neater
## Solving for the derivative:
odes:=solve(ode,{diff(a(t),t)});
ode1,ode2:=op(odes[1]),op(odes[2]);
dsolve(ode1); # implicitly given
## Finding equilibria:
RHS:=subs(a(t)=a, rhs(ode1));
E:=solve(RHS=0,a); #Equilibria: Same as solving sol[1] for a(t)
nops([E]); # 3 as expected, but not all are real
evalf(eval(E,{c=1,w=2,z=3})); # Example
e1:=%[1]; # The real one
eval(sol[1],{a(t)=e1,c=1,w=2,z=3}); # basically 0=0 as it should be
## 
## Numerical solution using the parameters option:
res:=dsolve({ode1,a(0)=1},numeric,parameters=[c,w,z]);
res(parameters=[1,2,3]); # Setting parameters of choice.
plots:-odeplot(res,[t,a(t)],-1..10,caption=typeset("Reaching equilibrium ",e1," in finite time"));
## The equilibrium is reached in finite time because of a sqrt(e1-a(t)) in ode1.
## We could find the time of arrival to e1 this way, if we want:
resE:=dsolve({ode1,a(0)=1},numeric,parameters=[c,w,z],events=[[sol[1],halt]]);
resE(parameters=[1,2,3]);
resE(8); # Time 5.9388290

Just name the inputs, as in
 

ode := 0.205*diff(u(t), t, t) + 0.53776*diff(u(t), t) + 7.49*u(t);
ics := u(0) = 0.6039948504, D(u)(0) = 0.2696963684;
dsolve([ode, ics]);

 

Begin by trying Y in this modified version where I only give the 2 procedures. Everything else is the same:
 

# Procedure to return values for specific time and parameter values.
u := proc( t, a, b ) 
	sol( 'parameters' = [ ':-a' = a, ':-b' = b ] ):
	try
		return eval( y( ':-t' ), sol( t ) ):  ### Changed
	catch:
		return 10000:
    end try:
end proc:


# Procedure to return objective function.
r := proc( a, b )
	local A, t, p, q, x, y:
	A := [ seq( u( t, a, b ), t=T ) ] -~ Y:  ### Changed
	return foldl( (p,q) -> p + q^2, 0, op(A) ):
end proc:

With X instead of Y Minimize cannot find an improved version.
A modified version allowing for several lists X,Y,Z, .. would consist in just replacing the zip part as follows:
 

`[]`~(X,Y) ## replacement for zip( (x,y) -> [x,y], X, Y )

This would work with no other changes in the original code for the 2 procedures. In the X,Y case the full line would be:
 

A := ListTools:-Flatten( [ seq( u( t, a, b ), t=T ) ] -~  `[]`~(X,Y) );

If only using one of them just remove the other, as in  `[]`~(Y).
One minor thing:
The return value in the catch error case ought to be of the same type as in the no error type:
i.e. [10000,10000] when using X,Y and [10000] when only using Y.

A note about the failure of Minimize when using X.
In the exact result for the system we see that x(t) only depends on the product a*b, whereas y(t) depends on both a and b independently.
This is going to be a problem for any numeric solver:
Try this version tailored to that situation:
 

ux := proc( t, ab ) #name ab. You could use any name of course
	sol( 'parameters' = [ ':-a' = ab, ':-b' = 1 ] ): # Arbitrarily setting b = 1.
	try
		return eval( [ x( ':-t' ) ], sol( t ) ):  
	catch:
		return [10000]:
    end try:
end proc:


# Procedure to return objective function.
rx := proc( ab ) #one input ab
	local A, t, p, q, x, y:
	A := ListTools:-Flatten( [ seq( ux( t, ab ), t=T ) ] -~  `[]`~(X) );  
	return foldl( (p,q) -> p + q^2, 0, op(A) ):
end proc:

rx(3);
# Find parameter values for bet fit.
Optimization:-Minimize( 'rx'(ab), ab=-1..1); # One variable ab

 

If you eliminate i(t) as done below your problem only depends on 3 parameters instead of 5.
That ought to help no matter what the method is.
 

restart;

Digits:=15:
eq1:=J*diff(theta(t),t,t)+b*diff(theta(t),t)=K*i(t);
eq2:=L*diff(i(t),t)+R*i(t)=V-K*diff(theta(t),t);

ICs:= theta(0)=0, D(theta)(0) = 0, i(0) = 0;
####
M:=ImportMatrix("F:/MapleDiverse/servo_open_ident_1_A.csv",source=csv):
M[1..10,..]; #Jump in V (col. 2) at start
M[-10..-1,..];
plot(M[..,1..2],labels=[t,V]); 
#### For V I shall use this:
V1:=Statistics:-NonlinearFit(a+b*sin(c*t+d),M[..,1..2],t,initialvalues={a=0,b=0.5,c=1/50,d=0});
plot(V1,t=0..50); p1:=%:
plot(M[..,1..2]); p2:=%:
plots:-display(p1,p2);
w1:=unapply(V1,t)~(M[..,1]):
w:=M[..,2]:
dw:=w-w1:
eval(V1,t=0);
dw[1];
(max,min)(dw[2..]); # Not bad
####
plot(M[..,[1,3]],size=[1000,default]); pd:=%:  # For later use kept in pd.
plot(M[..,[1,4]],size=[1000,default]); # Isn't used in the following
params:={J= 6.9347e-005, b = 2.2480e-004, K= 0.0094, R=0.2124, L=0.0015}; # Supplied by OP
####From eq1 and ICs follow that (D@@2)(theta)(0)=0. Will be used below.
####Eliminating i(t) by using integrating factor exp(R/L*t) used to solve first order linear odes:
S:=Diff(exp(R/L*t)*i(t),t)=exp(R/L*t)*rhs(eq2)/L;
##Operating on eq1:
map(Diff,exp(R/L*t)*convert(eq1,Diff)/K,t);
##Elimination of i(t):
subs(S,%);
ode:=value(%);
ode1:=L*J*op(solve(ode,{diff(theta(t),t,t,t)}));
ode2:=collect(ode1,diff,factor);
##We see that this third order ode depends on the following groups only: 
TR:={c1=(J*R+L*b)/(L*J), c2=(K^2+R*b)/(L*J), c3= K/(L*J)};
cparams:=eval(TR,params); #Op's params translated to c1,c2,c3
sbs:=solve(TR,{J,K,b});
ode3:=subs(sbs,ode2);
solve(ode3,{diff(theta(t),t,t,t)});
ODE:=eval(op(%),V=V1); # The 3rd order ode for theta(t).
##The initial conditions:
ics:=theta(0) = 0,D(theta)(0) = 0,(D@@2)(theta)(0)=0;
res:=dsolve({ODE,ics},numeric,parameters=[c1,c2,c3],output=listprocedure,maxfun=10^6);
th:=subs(res,theta(t));
## Q sets the parameters and computes the sum of squares of the residuals:
Q:=proc(c1,c2,c3) option remember; local j,ssq;
  if not [c1,c2,c3]::list(realcons) then return 'procname(_passed)' end if;
  try
    res(parameters=[c1,c2,c3]);
    ssq:=add( (th(M[j,1])-M[j,3])^2,j=1..5001)
  catch:
    ssq:=10^15
  end try;
  ssq
end proc:
##
t0:=time[real]():
RES:=Optimization:-NLPSolve(Q(c1,c2,c3),initialpoint=cparams);
time[real]()-t0; # 640s in my case.
##Finished with the error: No improved point could be found.
##But let us look at the remember table for Q:
T:=op(4,eval(Q)):
LL:=[entries(T,pairs)]:
LL:=remove(hastype,LL,{function,HFloat(undefined)}):
LS:=sort(LL,key=rhs):
LS[1];
LS[-1]; # A negative value for c1
res(parameters=([c1,c2,c3]=~[lhs(LS[1])]));
plots:-odeplot(res,[t,theta(t)],0..50,color=blue); pm:=%:
plots:-display(pd,pm,legend=[Data,Model]); #Not bad

The final plot of theta:

The c-values found were: [c1 = 0.07085248546343231, c2 = 1072.44657764141, c3 = 91326.7791941768]
(the result of [c1,c2,c3]=~[lhs(LS[1])]; )

 

Why not begin by solving this simple linear system exactly?
 

restart;
J:= 0.01;b:=0.1;K:= 0.01;R:=1;L:=1;
eq1:=J*diff(theta(t),t,t)+b*diff(theta(t),t)=K*i(t);
eq2:=L*diff(i(t),t)+R*i(t)=V-K*diff(theta(t),t);

ICs:= theta(0)=0, D(theta)(0) = 0, D(theta)(0) = 0, i(0) = 0:

res:=dsolve({eq1,eq2,ICs});

Then take it from there?

I suppose you might want to remove all equations of the form a = a ?
Then this would do it:
 

remove(x->(lhs-rhs)(x)=0,EQI__4);

or shorter:
 

remove( lhs-rhs = 0, EQI__4);

 

The result of splitting and a change of variable below seems alright, but the first result without splitting is I'm afraid a result of carelessness.
 

restart;
J := Int(cos(2*x)/(1+2*sin(3*x)^2), x = 0 .. Pi);
J1:=combine(J);
value(J1); # returns unevaluated (except with int instead of Int)
IntegrationTools:-Change(J1,cos(2*x)=z); # result right, but is it due to carelessness?
S:=IntegrationTools:-Split(J1,[Pi/4,Pi/2,3*Pi/4]);
L:=convert(S,list); # To see the steps
map(IntegrationTools:-Change,L,cos(2*x)=z);
convert(%,`+`);

 

Your system satisfies the requirements of uniqueness for an initial value problem, and that is what you got.
So since r[1] = V = q = 0 is a solution, it is the unique solution.

@siddikmaple Try Optimization:-LSSolve instead:
 

CSTR:={p[b], p[s], p[0, b], p[0, c], p[0, p], p[0, s], p[1, b], p[1, c], p[1, p], p[1, s], p[2, b], p[2, c], p[2, p], p[2, s], p[3, b], p[3, c], p[3, p], p[3, s], r[0], r[1], r[2], r[3], t[0], t[1], t[2], t[3]}=~(0..1);
##
residuals:=(lhs-rhs)~([e1,e2,e3,e4,e5,e6,e7,e8,e9,e10,e11,e12,e13,e14,e15,e16,e17,e18,e19,e20, e21,e22, e23, e24, e25, e26]);
##
ansLS:=Optimization:-LSSolve( residuals,op(CSTR),iterationlimit=10000);

The answer given is:

Continue with using that result as an initialpoint in SolveEquations:
 

DirectSearch:-SolveEquations(residuals,CSTR,initialpoint=ansLS[2]);

The answer is the same. Notice that LSSolve returns half of the sum of squares and gets:
0.563125933394723333e-1
whereas SolveEquations returns the sum of squares and finds:
.112625186678945
So there is basic agreement.

@siddikmaple You cannot use [ ] or { } instead of parentheses ( ) .
[ ] is used for lists and also for selection. { } is used for sets.
You have 4 cases , where [ ] is used instead of ( ), and 5 cases where { } is used instead of ( ):

sys:={e1, e10, e11, e12, e13, e14, e15, e16, e17, e18, e19, e2, e20, e21, e22, e23, e24, e25, e26, e3, e4, e5, e6, e7, e8, e9}:
nops(indets(sys,list)); # 4
nops(indets(sys,set)); # 5

This must be corrected before having any chance of getting a result out of fsolve.

Your system can be rewritten as a first order system of 4 dependent variables.
I have changed your code some and made extensive use of elementwise operations (using the tilde ~):
 

restart;
A := Matrix(2, 2, [0, -1, 1, -1]);
X := t-> <x(t), y(t)>;
g := X(t)/(t+1)^2;
alpha := -1/4;
X0 := <1, 1>;
B := exp(-2*alpha*(u-t))*Matrix(2, 2, [0, 1, 1, 0]);
eq := diff~(X(t), t) = A.X(t)+int~(B.X(u), u = 0 .. t)+g;
eqs:=expand~(Equate(lhs(eq),rhs(eq)));
S:=[y1(t)=int(exp(u/2)*y(u), u = 0 .. t),x1(t)=int(exp(u/2)*x(u), u = 0 .. t)];
sys1:=subs(rhs~(S)=~lhs~(S),eqs);
sys2:=diff~(S,t);
SYS:={sys1[],sys2[]};
## x(0)=1 and y(0)=1 are given; The last two follow from S.
ics:={x(0)=1,y(0)=1,x1(0)=0,y1(0)=0}; 
res:=dsolve(SYS union ics,numeric);
plots:-odeplot(res,[x(t),y(t)],-0.9..6);

The answer to your actual question about stability I shall leave for you to ponder about.
The plot produced in the code is:

If you look at the answer you get with c>0 you see that e.g. BesselK(0, s/c) appears. For this to make sense c must be different from zero. Thus besides c::real you would need c<>0.  But notice also that BesselK(0,x) is imaginary for x < 0. Thus the answer given for c>0 will probably not be valid for c<0.

That you get the same answer when using c>=0 as when using c>0 is puzzling, though.
But equalities and their negations ( <>) are often ignored by 'assume'.

Try solving for the derivatives and the problem will be obvious:
 

restart;
bvp:={ds*(diff(i(x), x)) = exp(eta(x)), s(x)^3*i(x)*b*r+(1-s(x))^3*a*l*(diff(s(x), x))/s(x)^1.5 = (1-s(x))^3, diff(eta(x), x) = dk*(i(x)-1)/s(x)^p, i(0) = 0, i(1) = 1, s(0) = 1};
sys:=select(has,bvp,diff);
factor(solve(sys,diff~({i,s,eta}(x),x)));

You get for diff(s(x),x) the following:
diff(s(x), x) = s(x)^(3/2)*(s(x)^3*i(x)*b*r+s(x)^3-3.*s(x)^2+3.*s(x)-1.)/(a*l*(s(x)-1.)^3)
As you see the denominator has the factor: (s(x)-1.)^3 which clearly is bad news if you want s(0) = 1.
The numerator is zero at x = 0, since i(0)=0, but it can be written as

s(x)^(3/2)*(s(x)^3*i(x)*b*r+(s(x)-1.)^3)/(a*l*(s(x)-1.)^3);

and again as
s(x)^(3/2)*(s(x)^3*i(x)*b*r/(s(x)-1.)^3 + 1)/al;

Clearly the expression s(x)^3*i(x)*b*r/(s(x)-1.)^3 is very problematic since for this to have no singularity at x = 0 i(x) would have to behave like i(x) = O(x^3) at x = 0.

 

I don't think that you are doing anything wrong. But there is no solution, which is also clear by inspection of the LinearSolve result.
We see that s[2], s[5], and s[6] must be positive. But then -3*s[2]-s[5]-2*s[6] is negative.
 

solve({x[1]>=0, x[2]>=0, x[3]>=0, x[4]>=0, x[5]>=0, x[6]>=0},[s[1],s[2],s[3], s[4], s[5],s[6]]);

returns [ ].
This could be a bug in LinearMultivariateSystem that it is unable to handle "no solution" in a graceful way.
Try this:
 

restart;
M:=Matrix([[3, 0, 2, 2, 1, 1], [-1, -1, 0, -1, 0, -1], [0, 3, 0, 1, 1, 2], [3, 0, 1, 2, 0, 1], [-1, -1, -1, -1, -1, -1]]);
vec:= Vector[column](5, [3, 1, 0, 0, -2]);
x:=LinearAlgebra:-LinearSolve(M,vec, free = s);
## Debugging SimpleAnd
debug(SolveTools:-Utilities:-SimpleAnd);
SolveTools:-Inequality:-LinearMultivariateSystem({x[1]>=0, x[2]>=0, x[3]>=0, x[4]>=0, x[5]>=0, x[6]>=0},[s[1],s[2],s[3], s[4], s[5],s[6]]);

Notice the green line saying

"enter Utilities:-SimpleAnd, args = -1 = 0, 0 <= -1+2*s[2]+s[5] "

The problem seems to be the false statement -1 = 0 having no parameter.

I have submitted a bug report (SCR).

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