C_R

2592 Reputation

19 Badges

5 years, 20 days

MaplePrimes Activity


These are answers submitted by C_R

I was expecting a more knowledgeable user than me to reply.

Anyway, I assume that the classification is wrong based on the following.

The method homogeneous is not appropriate for second order ode's because it solves first oder ode's under which it is also listed.

`dsolve/methods`[1];
  [quadrature, linear, Bernoulli, separable, inverse_linear, 

    homogeneous, Chini, lin_sym, exact, Abel, pot_sym]

Contrarily

dsolve(ode,[dAlembert])

solves the ode.

(Edit: there seems to be no homogenous,Class A method

`dsolve/methods`[1,high_degree];
dsolve(ode,y(x),%);
dsolve(ode,y(x),%%[-2..-1]);

Only two mehtods provide solutions)

Further, I do not see a match between


and

where I cannot convert the right hand side to an expression that exclusively depends on y(x)/x.

(Edit: It can be converted. See above reply to Carl.)

Also, the term homegeneous as introduced in the vocabulary in Advanced Engineering Mathematics with Maple

"A linear equation is homogeneous if, when the unknown function and all its derivatives are placed on the left side of the equation, the right side is then zero.  A linear equation is nonhomogeneous (inhomogeneous) otherwise.  Thus, the terms homogeneous and nonhomogeneous rightly apply to linear equations."

does not match the ode. Neither the right hand side is zero nor the left hand side represents a linear ode.

By the way, since there is no 1:1 correspondance between ode types and methods I used to look-up the methods

indices(`dsolve/methods`);
entries(`dsolve/methods`)

The code snippet is incomplete. I have higlighted in orange what I have added to understand what is going on.

I have stopped at the point where I found the name i not beeing assinged to a value.

 

NULL

with(geometry);
Sq := proc(n::integer)
local aS, oS, aC, oC, s, dr, pc, u;
aS := -i/n;
print(aS);
oS := sum(1/s, s = 1 .. n);
aC := 1/2*aS;
oC := oS - 1/2*1/((n + 1)*n);
point(S, aS, oS); point(C, aC, oC);
print(S);
MakeSquare(K, [S, 'center' = C]);
u := (x, i) -> sum(exp(-x*k)/k, k = 1 .. i);
pc := plot(u(x, n), x = 0 .. 4, color = green);
print([K]);
dr := draw([K]);
plots:-display({dr, pc});
end ;Sq(1);

[Apollonius, AreCollinear, AreConcurrent, AreConcyclic, AreConjugate, AreHarmonic, AreOrthogonal, AreParallel, ArePerpendicular, AreSimilar, AreTangent, CircleOfSimilitude, CrossProduct, CrossRatio, DefinedAs, Equation, EulerCircle, EulerLine, ExteriorAngle, ExternalBisector, FindAngle, GergonnePoint, GlideReflection, HorizontalCoord, HorizontalName, InteriorAngle, IsEquilateral, IsOnCircle, IsOnLine, IsRightTriangle, MajorAxis, MakeSquare, MinorAxis, NagelPoint, OnSegment, ParallelLine, PedalTriangle, PerpenBisector, PerpendicularLine, Polar, Pole, RadicalAxis, RadicalCenter, RegularPolygon, RegularStarPolygon, SensedMagnitude, SimsonLine, SpiralRotation, StretchReflection, StretchRotation, TangentLine, VerticalCoord, VerticalName, altitude, apothem, area, asymptotes, bisector, center, centroid, circle, circumcircle, conic, convexhull, coordinates, detail, diagonal, diameter, dilatation, directrix, distance, draw, dsegment, ellipse, excircle, expansion, foci, focus, form, homology, homothety, hyperbola, incircle, inradius, intersection, inversion, line, medial, median, method, midpoint, orthocenter, parabola, perimeter, point, powerpc, projection, radius, randpoint, reciprocation, reflection, rotation, segment, sides, similitude, slope, square, stretch, tangentpc, translation, triangle, vertex, vertices]

 

proc (n::integer) local aS, oS, aC, oC, s, dr, pc, u; aS := -i/n; print(aS); oS := sum(1/s, s = 1 .. n); aC := (1/2)*aS; oC := oS-(1/2)/((n+1)*n); geometry:-point(S, aS, oS); geometry:-point(C, aC, oC); print(S); geometry:-MakeSquare(K, [S, 'geometry:-center' = C]); u := proc (x, i) options operator, arrow; sum(exp(-x*k)/k, k = 1 .. i) end proc; pc := plot(u(x, n), x = 0 .. 4, color = green); print([K]); dr := geometry:-draw([K]); plots:-display({dr, pc}) end proc

 

-i

 

S

 

[K]

 

Error, (in geometry:-draw) non-numeric coordinate encountered, cannot determine plot view

 

NULL


 

Download incomplete_snippet.mw

Maple behaves deterministically until it applies the ParallelRisch method. I encountered a Maple crash while investigating this using the attached file. After further attempts on older files using printlevel commands, I experienced a severe system crash that necessitated two reboots to get my system, including mouse control, (believe it or not) back working.

Regarding your question "what is going on," I can offer a potential explanation:

Maple consistently enters ParallelRisch but appears to have challenges handling the integrand. Assuming "parallel" means Maple creates subtasks processed through multithreading (potentially on separate cores), variations in completion timings between subtasks could be a possibility. This might explain the observed behavior where timelimit stops processing in two instances, while the maple server completely freezes in another. I could test this hypothesis running Maple on a single core, which so far only produced the PDEtools/NumerDenom error and only one time a total freeze (this is not conclusive enough).

Other factors on CPU time variations might include additional tasks occasionally performed by mserver.exe alongside the integration process.

Since timelimit is a built-in routine, further investigation into it are not possible from my side.

I don't know whether it makes sense to invest more time in the combination of this special integrant and ParallelRisch. My understanding of ParallelRisch is that is does not more than Risch does. If so, ParalleRisch should have returned the input unevaluated as Risch does.

Can the below be reproduced? And OS different form Windows would be of interest. In any case, better to save your work before trying...

infolevel[int]:=5;

5

(1)

expr:=-4*(1-exp(I*x))^(-4*x)*(exp(I*x)+1)^(4*x)*exp(4*I*(polylog(2,exp(I*x)))-polylog(2,-exp(I*x) ))*csc(x)*x*(tan(x)^2-1);

-4*(1-exp(I*x))^(-4*x)*(exp(I*x)+1)^(4*x)*exp((4*I)*polylog(2, exp(I*x))-polylog(2, -exp(I*x)))*csc(x)*x*(tan(x)^2-1)

(2)

int(expr,x,method=ParallelRisch);
(* This command makes mserver "unresponsive" after a while: interrupt from the GUI does not work, Time in the status bar is not updated, private bytes = system alloacted memory jump 64 Gb, physical memory (visible in Windows task manager) still piling up*)

int: Beginning integration with _EnvContinuous=_EnvContinuous, _EnvAllSolutions=_EnvAllSolutions, and _EnvCauchyPrincipalValue=_EnvCauchyPrincipalValue.

 

Degree to be used when constructing the parallel integration system = 11

 

(* this input froze Maple GUI*, output was backed up before excecution of input
- might required sereval attempts to reproduce, which was possible*)
printlevel:=10;
st:=time[real]():
timelimit(60,int(expr,x,method=ParallelRisch));
print("time taken ",time[real]()-st);

int: Beginning integration with _EnvContinuous=_EnvContinuous, _EnvAllSolutions=_EnvAllSolutions, and _EnvCauchyPrincipalValue=_EnvCauchyPrincipalValue.

 

Degree to be used when constructing the parallel integration system = 11

 

Error, (in PDEtools/NumerDenom) time expired

 

"time taken ", 82.701

(3)

NULL

Download untitled2150_MAS.bak.mw

Its not fully clear what the unknows are that you want dsolve to solve for.

nops(sys)
                               12

tells us that you have 12 equations.

indets(sys);
 /                                      (1/2)   d           
{ exp, omicron, t, theta, (42.3 - theta)     , --- i[h](t), 
 \                                              dt          

   d            d            d            d           
  --- i[v](t), --- j[v](t), --- r[h](t), --- s[h](t), 
   dt           dt           dt           dt          

   d           
  --- s[v](t), 
   dt          

     /                                                  2\  
  exp\0.02366863905 (-0.0464 omicron + 1.046 theta - 23) /, 

     /                                                     4\  
  exp\0.00002055229087 (0.01 omicron + 1.01 theta - 21.211) /, 

  exp(0.0054 theta + 0.6737), i[h](t), i[v](t), j[v](t), r[h](t), 

                  \ 
  s[h](t), s[v](t) }
                  / 


lists additional to 12 unknown functions of t the parameters omicron and theta. These are not functions of t as you have defined them in the first two input lines. As such sys cannot be solved numerically.

Also: The frist two input lines do not use the assignement operator to define theta and omicron.

This works as well

map(X->(odetest(X,ode) assuming x>0),[sol])

The parenthesis group terms. In that way the solutions are mapped onto an expression that inculdes the assumption.

Afterwards odetest is excecuted with the assumptions.

Without parentheses the parser does not spot the second argument.

This one works as well

eval(map(X->'odetest(X,ode) assuming x>0',[sol]));

The reason "why Maple gives this internal error" is that Maple applies the procedure SolveSeparable which provides a result that gives an error because it computes ln(0) for the IC y(0)=1.

Only for this IC the result does not work.

In a similar case the Student,ODEs package had an appropriate method on board but the case (method not appropriate for a particular IC) could be solved by fixing the false classification "exact" to "NONE".
For this case it seems that there is no Student,ODEs method.

dsolve solves this IVP differently with an implicit method.

infolevel[dsolve]:=5;
dsolve({ode,ic},implicit);
solve(%,{y(x)})[];
dsolve({ode,ic});

Since I don't see that an implicit method is implemented in ODESteps, I'm curious if it's really a bug or a current limitation of ODESteps.

You can extract the polygons (vertical bars) from the plot structure which decribe the xy data.
The polygons have four corners. You probably want not xy data from all corners.
In the below I have extracted the very left x coordinate of the polygons and added the most right corner of the histogram.

[op(op(1, Q))];
%[1 .. nops(%) - 1];
seq(%[i][1][1], i = 1 .. nops(%)), %[nops(%)][2][1];

You can do the same for y coordinates or pairs of x,y whatever you are interesed in.

This is a crude example how you could filter subexpressions still using indets

indets(eqn) minus indets(eqn,name):
subexpr:=map(op,%) minus indets(%,fraction):# remove of root terms - instead of filtering for `^`
rpl:=[A,B,C,`&D;`,E,F]:
Subs2:=[seq(subexpr[i]=rpl[i], i=1..nops(subexpr))]:
simplify(subs(Subs2,eqn));

This could be refined and put into a procedure.( I am not quite sure if you are looking for something like that.)

ODESteps does not apply an adequate method although it has one on board:
Student:-ODEs:-Solve:-FirstOrderLinear. Instead "homotopy method" is applied.

When tracing ODEstep I got the impression that ODEstep tries to follow dsolve closely and adds comments step by step. Programming this is probably non trivial because essential steps in the execution of a sophisticated command (that was initially not designed for the purpose of education) have to be identified and attributed to a mathematical correct instruction to a student.

Assuming no bug in the program flow, the answer to your question “Why ODEStep is unable” could be: FirstOrderLinear is not implemented yet. This is likely because it is not listed here under “Step-by-Step Computation”.
I guess we will see it implemented in the future.

 

Update:
The old link did not work because of a special character in the file name.
first_order_quadrature_ode_with_odesteps.mw

(now v1749 can compute the integral using two methods )

Maple distinguishes between many units because they either are in use or they have been. For some units the context is important. ?Units,annotations provides some examples. If you find them strange or exotic it’s probably because you use the SI system and do not have to deal with, for example, historic units.
Maple also provides energy conversions ?Units,EnergyConversions which seem to be useful in certain areas but are not known at all in others domains.

The units package has grown over the years and includes allot of options and functions users have requested but other users do not need. For this reason the package has become powerfull but not straighforward to use (at least it was for me).
Up to you to use Grays or not.

 

Personally, I use the default unit settings in combination with convert and simplify. The context panel is usefull too.
If you have any particular calculation or conversion problems consider uploading a file with


 

2*Unit('gray')

2*Units:-Unit(Gy)

(1)

Using the context panel and conversion to E0/em and erg

2*Units:-Unit(Gy)

2*Units:-Unit(Gy)

(2)

2*Units:-Unit(Gy)

2*Units:-Unit(Gy)

(3)

Now: Remove the context

convert(2*Units:-Unit(Gy), units, J/kg)

2*Units:-Unit(J/kg)

(4)

` `After that, the context pannel does not offer conversion to Gray anymore (but to other units)

2*Units:-Unit(J/kg)

2*Units:-Unit(J/kg)

(5)

NULL

Lets do a calculation

2*Units:-Unit(Gy)-2*Units:-Unit(J/kg)

2*Units:-Unit(Gy)-2*Units:-Unit(J/kg)

(6)

simplify(2*Units:-Unit(Gy)-2*Units:-Unit(J/kg))

0

(7)

NULL


 

Download Unit_conversion.mw

The problem with your algorithm is the break statement that terminates the inner loop without defining a condition for t[j+1].

When the outer loop continues and increments j from 2 to 3 then this statement

assigns t[3] to u2[0]. t[3] however has no numeric value assigned to and Maple starts evaluating symbolic expressions. This blows up the whole calculation without giving any usefull results.

 

Attached is your code with debugging statements that allow you to investigate line by line the computation.

?worksheet,interactive,debugger explains how to use it.

NULL

 

restart;

 

dbg:=proc ()

DBG> quit

showstat(dbg)


dbg := proc()
local HAM, TOL, x, N, h, M, t, alpha, f1, f2, f3, g1, g2, j, u1, u2, u3, v1, v2, i,
kf1, kf2, kf3, kg1, kg2, U;
   1   HAM := [1];
   2   Digits := 30;
   3   TOL := 1/1000;
   4   x[0] := 0;
   5   N := 5000;
   6   h := .1e-2;
   7   M := 20;
   8   t[0] := .6;
   9   alpha := 0.;
  10   f1 := (x, u1, u2, u3, v1, v2) -> u2;
  11   f2 := (x, u1, u2, u3, v1, v2) -> u3;
  12   f3 := (x, u1, u2, u3, v1, v2) -> u2^2-2*u1*u3-v1^2;
  13   g1 := (x, u1, u2, v1, v2) -> v2;
  14   g2 := (x, u1, u2, v1, v2) -> 2*u2*v1-2*u1*v2;
  15   for j from 0 while j <= M do
  16       u1[0] := 0;
  17       u2[0] := t[j];
  18       u3[0] := 0;
  19       v1[0] := .8;
  20       v2[0] := 0;
  21       for i from 0 while i <= N-1 do
  22           DEBUG(is((i = 221)*`and`(j = 2)));
  23           kf1[i,1] := f1(x[i],u1[i],u2[i],u3[i],v1[i],v2[i]);
  24           kf2[i,1] := f2(x[i],u1[i],u2[i],u3[i],v1[i],v2[i]);
  25           kf3[i,1] := f3(x[i],u1[i],u2[i],u3[i],v1[i],v2[i]);
  26           kg1[i,1] := g1(x[i],u1[i],u2[i],u3[i],v1[i],v2[i]);
  27           kg2[i,1] := g2(x[i],u1[i],u2[i],u3[i],v1[i],v2[i]);
  28           kf1[i,2] := f1(x[i]+1/2*h,u1[i]+1/2*h*kf1[i,1],u2[i]+1/2*h*kf2[
                 i,1],u3[i]+1/2*h*kf3[i,1],v1[i]+1/2*h*kg1[i,1],v2[i]+1/2*h*
                 kg2[i,1]);
  29           kf2[i,2] := f2(x[i]+1/2*h,u1[i]+1/2*h*kf1[i,1],u2[i]+1/2*h*kf2[
                 i,1],u3[i]+1/2*h*kf3[i,1],v1[i]+1/2*h*kg1[i,1],v2[i]+1/2*h*
                 kg2[i,1]);
  30           kf3[i,2] := f3(x[i]+1/2*h,u1[i]+1/2*h*kf1[i,1],u2[i]+1/2*h*kf2[
                 i,1],u3[i]+1/2*h*kf3[i,1],v1[i]+1/2*h*kg1[i,1],v2[i]+1/2*h*
                 kg2[i,1]);
  31           kg1[i,2] := g1(x[i]+1/2*h,u1[i]+1/2*h*kf1[i,1],u2[i]+1/2*h*kf2[
                 i,1],u3[i]+1/2*h*kf3[i,1],v1[i]+1/2*h*kg1[i,1],v2[i]+1/2*h*
                 kg2[i,1]);
  32           kg2[i,2] := g2(x[i]+1/2*h,u1[i]+1/2*h*kf1[i,1],u2[i]+1/2*h*kf2[
                 i,1],u3[i]+1/2*h*kf3[i,1],v1[i]+1/2*h*kg1[i,1],v2[i]+1/2*h*
                 kg2[i,1]);
  33           kf1[i,3] := f1(x[i]+1/2*h,u1[i]+1/2*h*kf1[i,2],u2[i]+1/2*h*kf2[
                 i,2],u3[i]+1/2*h*kf3[i,2],v1[i]+1/2*h*kg1[i,2],v2[i]+1/2*h*
                 kg2[i,2]);
  34           kf2[i,3] := f2(x[i]+1/2*h,u1[i]+1/2*h*kf1[i,2],u2[i]+1/2*h*kf2[
                 i,2],u3[i]+1/2*h*kf3[i,2],v1[i]+1/2*h*kg1[i,2],v2[i]+1/2*h*
                 kg2[i,2]);
  35           kf3[i,3] := f3(x[i]+1/2*h,u1[i]+1/2*h*kf1[i,2],u2[i]+1/2*h*kf2[
                 i,2],u3[i]+1/2*h*kf3[i,2],v1[i]+1/2*h*kg1[i,2],v2[i]+1/2*h*
                 kg2[i,2]);
  36           kg1[i,3] := g1(x[i]+1/2*h,u1[i]+1/2*h*kf1[i,2],u2[i]+1/2*h*kf2[
                 i,2],u3[i]+1/2*h*kf3[i,2],v1[i]+1/2*h*kg1[i,2],v2[i]+1/2*h*
                 kg2[i,2]);
  37           kg2[i,3] := g2(x[i]+1/2*h,u1[i]+1/2*h*kf1[i,2],u2[i]+1/2*h*kf2[
                 i,2],u3[i]+1/2*h*kf3[i,2],v1[i]+1/2*h*kg1[i,2],v2[i]+1/2*h*
                 kg2[i,2]);
  38           kf1[i,4] := f1(x[i]+h,h*kf1[i,3]+u1[i],h*kf2[i,3]+u2[i],h*kf3[i
                 ,3]+u3[i],h*kg1[i,3]+v1[i],h*kg2[i,3]+v2[i]);
  39           kf2[i,4] := f2(x[i]+h,h*kf1[i,3]+u1[i],h*kf2[i,3]+u2[i],h*kf3[i
                 ,3]+u3[i],h*kg1[i,3]+v1[i],h*kg2[i,3]+v2[i]);
  40           kf3[i,4] := f3(x[i]+h,h*kf1[i,3]+u1[i],h*kf2[i,3]+u2[i],h*kf3[i
                 ,3]+u3[i],h*kg1[i,3]+v1[i],h*kg2[i,3]+v2[i]);
  41           kg1[i,4] := g1(x[i]+h,h*kf1[i,3]+u1[i],h*kf2[i,3]+u2[i],h*kf3[i
                 ,3]+u3[i],h*kg1[i,3]+v1[i],h*kg2[i,3]+v2[i]);
  42           kg2[i,4] := g2(x[i]+h,h*kf1[i,3]+u1[i],h*kf2[i,3]+u2[i],h*kf3[i
                 ,3]+u3[i],h*kg1[i,3]+v1[i],h*kg2[i,3]+v2[i]);
  43           u1[i+1] := u1[i]+1/6*h*(kf1[i,1]+2*kf1[i,2]+2*kf1[i,3]+kf1[i,4]
                 );
  44           u2[i+1] := u2[i]+1/6*h*(kf2[i,1]+2*kf2[i,2]+2*kf2[i,3]+kf2[i,4]
                 );
  45           u3[i+1] := u3[i]+1/6*h*(kf3[i,1]+2*kf3[i,2]+2*kf3[i,3]+kf3[i,4]
                 );
  46           v1[i+1] := v1[i]+1/6*h*(kg1[i,1]+2*kg1[i,2]+2*kg1[i,3]+kg1[i,4]
                 );
  47           v2[i+1] := v2[i]+1/6*h*(kg2[i,1]+2*kg2[i,2]+2*kg2[i,3]+kg2[i,4]
                 );
  48           if abs(v1[i+1]-1) <= TOL then
  49               N := i+1;
  50               print("Value of v = ",[(i+1)*h, v1[i+1], u2[N]]);
  51               if abs(u2[i+1]) <= TOL then
  52                   print("Value of f' = ",u2[i+1]);
  53                   break
                   else
  54                   if j = 0 then
  55                       t[j+1] := t[j]+(alpha-u1[i+1])/(i+1)/h;
  56                       U[j] := u2[i+1];
  57                       print("Value of u1[j] = ",%)
                       else
  58                       U[j] := u2[i+1];
  59                       print("Value of u1[j] = ",%);
  60                       t[j+1] := t[j]-(U[j]-alpha)*(t[j]-t[j-1])/(U[j]-U[j
                             -1]);
  61                       print("Value of t = ",%)
                       end if
                   end if
               end if
           end do;
  62       print(j,i)
       end do
end proc
 

 

stopat(dbg,49);

[dbg]

(1)

dbg()

2, 222

(2)

NULL


 

Download debug.mw

The reason is this line of code that (I assume) determines the integration constant

Since the solution contains Pi and c__1 as a name (one name too much)

solve throws and error