13613 Reputation

19 years, 253 days

Piecewise with method=classical...

@mmcdara Classical methods in dsolve (like method=classical[rk4]) use fixed stepsize and has no error control. That allows us to use your original piecewise version mentioned in your original question.
Results appear good with stepsize = 1e-4.

```# Using P_DATA100 the system is:
sys100 := {10*diff(v(t), t) = 900*t - 1000 - 10000*x(t) - 100*piecewise(0 < v(t), 1, v(t) < 0, -1, 0), diff(x(t), t) = v(t), v(10/9) = 0, x(10/9) = 0};
#Now using Runge-Kutta 4:
res:=dsolve(sys100,numeric,method=classical[rk4],stepsize=0.0001);
odeplot(res,[t,x(t)],T__0..2);```

A successful events version...

@mmcdara I have attached a successful events version.
I will say, however, that the classical approach is simpler and real fast with stiff=true (i.e. rosenbrock).
The discrete variable f(t) used here is of type float and will take on 3 values: -1, 0, or 1.
Events are triggered by v(t) +eps =0, v(t)-eps = 0, and v(t)=0.
For both methods stiff=true is chosen although it probably has no effect in the events case.

 > restart;
 > with(plots):
 > # x(t)        body displacement # v(t)        body velocity # M__p     body mass # V__p     body volume # P__0      pre-stress of the spring # K           spring stiffness # C           solid-solid friction coefficient # rho__h   density of the fluid the body is immersed in # # gamma__longi(t) longitudinal acceleration gamma__longi := t -> A*t;
 > ####################################
 > # The classical approach:
 > equationsC := {                diff(x(t), t) = v(t),                M__p * diff(v(t), t) = (M__p - V__p*rho__h) * gamma__longi(t)                                       - (P__0 + K*x(t))                                       - C*tanh(v(t)/v__0)             };
 > EVENTS:=[[v(t)+eps,f(t)=-1],[v(t)-eps,f(t)=1],[v(t),f(t)=0]];
 > equationsE := {                diff(x(t), t) = v(t),                M__p * diff(v(t), t) = (M__p - V__p*rho__h) * gamma__longi(t)                                       - (P__0 + K*x(t))                                       -f(t)*C             };
 > # take-off time # (obtained by equating to 0 the rhs of the 2nd edo after having set x(t)=v(t)=0 and removed the friction term) rest_state := eval(equationsC[1], [x(t)=0, v(t)=0]); T__off     := solve(eval(rest_state, C=0), t);
 > # initial conditions IC_C := { x(T__off) = 0, v(T__off) = 0}; IC_E := { x(T__off) = 0, v(T__off) = 0, f(T__off) = 0}; # full systems sysC := equationsC union IC_C; sysE := equationsE union IC_E;
 > eps:=1e-4: v__0:=1e-7:
 > SOL_C:=dsolve(sysC,numeric,parameters=[A, M__p, V__p, rho__h, P__0, K, C],stiff=true); ## stiff SOL_E:=dsolve(sysE,numeric,parameters=[A, M__p, V__p, rho__h, P__0, K, C],                    discrete_variables=[f(t)::float],                    events=EVENTS,maxfun=10^6, stiff=true  #stiff doesn't help much here                 );
 > P_DATA0 := [             A      = 100,             M__p   = 10,             V__p   = 1,             rho__h = 1,             P__0   = 1000,             K      = 10000,             C      = 0           ]: T__0 := eval(T__off, P_DATA0);
 > SOL_C(parameters=P_DATA0): Frictionless_C := odeplot(SOL_C, [t, x(t)], t=T__0..2, color=blue);
 > SOL_E(parameters=P_DATA0): Frictionless_E := odeplot(SOL_E, [t, x(t)], t=T__0..2, color=blue);
 > P_DATA100 := [             A      = 100,             M__p   = 10,             V__p   = 1,             rho__h = 1,             P__0   = 1000,             K      = 10000,             C      = 100           ]: T__0 := eval(T__off, P_DATA100);
 > SOL_C(parameters=P_DATA100): t0:=time(): StrongFriction_C := odeplot(SOL_C, [t, x(t)], t=T__0..2, color=red): time()-t0; # 0.015
 > display(Frictionless_C, StrongFriction_C); pC:=%:
 > SOL_E(parameters=P_DATA100): t0:=time(): StrongFriction_E := odeplot(SOL_E, [t, x(t)], t=T__0..2, color=red): time()-t0; # 2.250s
 > display(Frictionless_E, StrongFriction_E); pE:=%:
 > display(pC,pE);
 >

pC and pE look the same and overlap totally in the last display with eps=1e-4.
With eps = 1e-3 computation time is much smaller and overlap is still good.
With eps = 1e-2 computation time is even smaller, but there is clear visual disagreement between pC and pE.
odeplot(SOL_E,[t,f(t)], T__0..2,numpoints=10000) will show you that an enormous amount of events takes place, so no wonder it takes time.

Time...

@Muhammad Usman Commenting out what may interfere with the timing (such as remembering values from the plots) and omitting printing the loop I get 0.187s.
fnum can be seen afterwards simply by eval(fnum);

```restart;

U := q^6*sin(p);
alpha := 1/3;
M1 := 3;
M2 := 3;
A:=Int(Int(U/((x^rho - p^rho)^alpha*(y^rho - q^rho)^alpha), p = 0 .. x), q = 0 .. y);

Q:=proc(r,xx,yy) evalf(eval(A,{'rho'=r,x=xx,y=yy})) end proc;
##plot3d(Q(1,x,y),x=0..1,y=0..1); # Test
## An animation in rho with only 6 frames
##plots:-animate(plot3d,[Q(rho,x,y),x=0..1,y=0..1],rho=1..6,frames=6);

## The revised loop:
t0:=time():
##printlevel:=3:
for i to 6 do
for i11 to M1 do
for j11 to M2 do
fnum[i, i11, j11] := Q(i,i11/M1,j11/M2)
end do
end do
end do:
time()-t0;
# time = 0.187s
eval(fnum);
```

Another version...

@mmcdara I just tried these events:

```events=[[[0,v(t)>0],b(t)=1],[[0,v(t)<0],b(t)=0]]
```

They seem to produce exactly what you have in your event version.

Would like to...

@mmcdara I would like to, but haven't seen a way out yet.

ConvertIn...

When I do ?ConvertIn the help page for modp1 opens.
On that page it says:
"The modp1 function accepts the following functions, whose names are also protected global names."

Then follows 5 colums of names some of which are in black (i.e. no link to help page). Among those is ConvertIn.
Its meaning is understood by modp1 and is described on that help page. In fact ConvertIn may just be nothing but a protected name used in a call to modp1 to tell it what to do.
Some of those with an actual link (the blue ones) like Add just brings you to the help page for add. That must be a mistake.
How Add is used is beyond me.

The usual way...

A common way of solving this kind of ode is to multiply by diff(y(x),x) on both sides and integrate like this:

```restart;
ode:=2*diff(y(x),x\$2)=exp(y(x));
IC:=y(0)=0,D(y)(0)=1;
diff(y(x),x)*ode;
ode1:=map(int,%,x)+(0=C1);
eval[recurse](convert(ode1,D),{x=0,IC});
## So C1 = 0.
ode2:=eval(ode1,C1=0);
odes:=solve(ode2,{diff(y(x),x)});
ode3:=op(odes[1]);
sol1:=dsolve(ode3);
eval[recurse](sol1,{x=0,IC});
## so _C1 = -2.
sol:=eval(sol1,_C1=-2);
```

Is this the ode?...

Which one of these two (if any of them) did you mean:

```ode1:=diff(U(x), x\$2) + ((z+ I*y)^2/k1^2 -k2^2 + ((z+I*y)/k3)*sech(x)^2)* U(x);
ode2:=diff(U(x), x\$2) + ((z+ I*y)^2/(k1^2 -k2^2) + ((z+I*y)/k3)*sech(x)^2)* U(x);
```

Maple 2020 solves both symbolically.

A procedure and a revised loop:...

1. You can make a very simple procedure (Q below) to return values of the integral for given x and y values.
2. In your original triple loop you didn't keep all values computed. fnum[i] was redefined for every new values of i11 and j11.

```restart;

U := q^6*sin(p);
alpha := 1/3;
M1 := 3;
M2 := 3;
A:=Int(Int(U/((x^rho - p^rho)^alpha*(y^rho - q^rho)^alpha), p = 0 .. x), q = 0 .. y);
Q:=proc(r,xx,yy) evalf(eval(A,{'rho'=r,x=xx,y=yy})) end proc;
plot3d(Q(1,x,y),x=0..1,y=0..1); # Test
## An animation in rho with only 6 frames
plots:-animate(plot3d,[Q(rho,x,y),x=0..1,y=0..1],rho=1..6,frames=6);
## The revised loop:
printlevel:=3:
for i to 6 do
#rho := i; #Not used
for i11 to M1 do
for j11 to M2 do
fnum[i, i11, j11] := Q(i,i11/M1,j11/M2)
end do
end do
end do;
```

What exactly did you do?...

@Muhammad Usman Could you please give the full code of what you tried to do. Perhaps upload a worksheet?

Numerical integration...

@Muhammad Usman Numerical integration requires x and y to be given numerical values.

So here is an example where x = 1 and y = 2:

```for i from 1 to 6 do
rho := i;
fnum[i] := evalf(eval(Int(Int(U/((x^rho-p^rho)^alpha*(y^rho-q^rho)^alpha), p = 0 .. x), q = 0 .. y),{x=1,y=2}));
end do;
```

Compare that with:

`evalf(eval([seq](f[i],i=1..6),{x=1,y=2}));`

Revised version...

@Muhammad Usman I have revised my code above so U now is the revised U.
If you don't reverse the order of integration f[4] is not completed and other results are messier.

showstat...

@Anthrazit By looking at the codes for the Standard and Simple versions of min (and by debugging both):

```showstat(Units:-Standard:-min);
showstat(Units:-Simple:-min);
debug(Units:-Standard:-min);
debug(Units:-Simple:-min);

```

you will see that Units:-Simple:-min among several return statements has a return statement at the bottom

return min(op(arguments))
where the local 'arguments' is just the list of arguments passed to Units:-Simple:-min.
What :-min gets as arguments is easily seen by doing debug(:-min);
So after a restart you can try:

```restart;
min(15*Units:-Unit(m), 0);  # 0
max(15*Units:-Unit(m), 0); # 15*Units:-Unit(m)
```

Parameters...

@srikanthan It might help us if we knew the values of the parameters excepting of course R and thetaw, i.e.
{P, Pr, Q, Qc, S0, S1, S2, S4, c, lambda, n}.
Or at least an example of reasonable values.
Then we could make a procedure taking the value (Dth0) of D(theta)(0) and thetaw as input and R as output using as an extra boundary condition D(theta)(0) = Dth0.

n = 1...

If you insist on the initial conditions in dsys6 then it follows that n = 1. Simply use the connection between x and y given in the first equation together with the initial conditions.

I suppose then that there is some mistake here?

 First 22 23 24 25 26 27 28 Last Page 24 of 229
﻿