## 13222 Reputation

12 years, 354 days

## Rule number 1...

Use the big green up-arrow in the Mapleprimes toolbar

## You refuse...

The ode you define in your worksheets as  'de1' only contains the dependent variable 'f(x)'. The ode you define in your worksheets as  'de2' contains the dependent variables 'f(x)' and 'g(x)'.

When you process the ode 'de1', to produce 'DE1', the latter contains the dependent variables

{b[0](x), b[1](x), b[2](x), b[3](x), b[4](x), b[5](x), b[6](x), b[7](x), b[8](x), b[9](x), b[10](x)}

However when you process the ode 'de2', to produce 'DE2', the latter contains the dependent variables (check the one highlighted in red)

{f(x), c[0](x), c[1](x), c[2](x), c[3](x), c[4](x), c[5](x), c[6](x), c[7](x), c[8](x), c[9](x), c[10](x)}

The mere existence of the undetermined function 'f(x)' in the definition of 'DE2' is the source of your problem.

See the attached for proof - it makes no attempt to solve anything - it just lists the indeterminates. Until you fix the existence of the indeterminate 'f(x)' in the definition of 'DE2' you are going to get precisely nowhere

 > restart:
 > Pr:=1:
 > de1 := (1 - p)*diff(f(x), x \$ 3) + p*(diff(f(x), x \$ 3) + 1/2*f(x)*diff(f(x), x \$ 2)):
 > de2 := (1 - p)*diff(g(x), x \$ 2)/Pr + p*(diff(g(x), x \$ 2)/Pr + 1/2*f(x)*diff(g(x), x)):
 > ibvc := f(0), D(f)(0), D(f)(5) - 1, g(0) - 1, g(5):   n := 10:   F := unapply(add(b[k](x)*p^k, k = 0 .. n), x):   G := unapply(add(c[k](x)*p^k, k = 0 .. n), x):
 > DE1 := series(eval(de1, f = F), p = 0, n + 1):   indets(DE1, function(name));
 (1)
 > DE2 := series(eval(de2, g = G), p = 0, n + 1):   indets(DE2, function(name));
 (2)
 >

## Have you tried...

actually looking at the system you are trying to solve?

The final execution group in your latest worksheet - ie this piece of code

```for k from 0 to n do
IBVC1 := select(has, CT, c[k]);
{coeff(DE2, p, k), op(IBVC1)};
slv1 := dsolve({coeff(DE2, p, k), op(IBVC1)});
c[k] := unapply(rhs(slv1), x);
end do```

For k=0, your code generates the ODE system

{c[0](0) - 1, diff(c[0](x), x, x), c[0](5)}

a second order ODE with two boundary conditions - which is fine, and Maple generates the solution

slv1 := c[0](x) = -x/5 + 1

For k=1, your code generates the ODE system

{diff(c[1](x), x, x) - f(x)/10, c[1](0), c[1](5)}

a second-order ODE with two boundary conditions, and a completely unknown function of the same independent variable f(x). What solution do you expect?

Maple returns a "formal" solution containing double integrals of the unknown function f(x)

## Not sure what point...

you are trying to make.

1. Your original code generates an invalid ODE system, which means that dsolve() will produce an error.
2. This update to your code produces a valid ODE system, which means that dsolve() has a reasobale chance of coming up with a solution.

These two cases, both for loop index k=0 are illustrated in the attached. One errors, for reasons I explained previously and one doesn't.

 > # # OP's original code asks for a solution of #   dsolve({1, diff(c[0](x), x, x)}); # # which will produce an error #
 > # # OP's new code asks for a solution of #   dsolve( {D(b[0])(5) - 1, diff(b[0](x), x, x, x), b[0](0), D(b[0])(0)});
 (1)
 >

## Use the ignore=true option...

as in the attached "toy" example

 > restart;   with(Statistics):   V:=Vector[column](9, [124.0, 130.0, 130.0, 119.0, 136.0, 118.0, undefined, 130.0, 95.0]);   Mean(V, ignore=true);
 (1)
 >

## I assume...

that you were trying to achieve what is shown in the attached? I don't think this approach can be (easily) extended beyond first differences.

 > restart;   a := 1: b := 5: h := 1: f := 1/x: N := (b-a)/h:   for i from 0 while i <= N do       x[i] := h*i+a;       y[i] := eval(f, x = x[i])   end do:   printf("_______________________________________________________________________________\n\n\n");     for i from 0 by 1 while i<=2*N do:         if   irem(i,2)=0       then printf("%2d%16.7f%16.7f\n",i,x[i/2],y[i/2]);       else printf("%2d%48.7f\n",i, y[(i+1)/2]-y[(i-1)/2]);       fi;   od;   printf("_______________________________________________________________________________\n")
 _______________________________________________________________________________  0       1.0000000       1.0000000  1                                      -0.5000000  2       2.0000000       0.5000000  3                                      -0.1666667  4       3.0000000       0.3333333  5                                      -0.0833333  6       4.0000000       0.2500000  7                                      -0.0500000  8       5.0000000       0.2000000 _______________________________________________________________________________
 >

@ check my original comments on the typos in your equation and fix them properly before doing anything else.

I'd be surprised if an analytic solution can be found for the resulting ODE, so you may have to be satisfied with a numeric one - which means that you will need values for all parameters and a couple of initial/bondary conditions:-(

## With a trivial modification...

the attached will show which months in a specified year have a Friday 13th. Output is in the form

[year, [list of months with Friday 13th]]

From which is pretty obvious that in 2022, the only month with a Friday 13th is month 5 - ie May!

 > restart;   with(Calendar):   fri:= yr-> local j:seq(`if`(DayOfWeek(yr, j, 13)=6,j, NULL), j=1..12):
 > # # So which months in the supplied year have a Friday 13th #   seq( [j, [fri(j)]], j=2000..2022);
 (1)
 >
 >

## You mean...

when you type

in Maple Flow - it doesn't work?

Because if this is true then all you have demonstrated ist that Maple Flow is unable to execute basic Maple commands - which is somewhat scary!

## Idle curiosity...

If you want to restrict the range variable to integers, then not only do you need the 'sample' option, but ialso to set the adaptive=false option. If the intent is to restrict the range variable to integers, why have a 'floor(n)' command in the function definition? Isn't this a bit - well - superfluous?

## The point I was trying to make...

was that in a plot() command, the plotting variable is continuous - it does not assume integer values. Kitonum circumvented this problem by plotting only points. This has the drawback that if you draw lines between the points (say by using style=pointline), then you will not get "vertical" lines, since x-values will always differ by one

I gave a (quick+dirty) workaround for this - and you can see the effect of superimposing my plot and Kitonum's in the attached. This gives an apparent "right shift to the curve on my plot. If this is a problem it is relatively easy to fix, see the final figure in the attached

 > restart;   f:=n->ceil(sqrt(4*floor(n)))-floor(sqrt(2*floor(n)))-1:   plot(f, 10..100, gridlines=false);
 > Points:=[[10, 2], [11, 2], [12, 2], [13, 2], [14, 2], [15, 2], [16, 2], [17, 3], [18, 2], [19, 2], [20, 2], [21, 3], [22, 3], [23, 3], [24, 3], [25, 2], [26, 3], [27, 3], [28, 3], [29, 3], [30, 3], [31, 4], [32, 3], [33, 3], [34, 3], [35, 3], [36, 3], [37, 4], [38, 4], [39, 4], [40, 4], [41, 3], [42, 3], [43, 4], [44, 4], [45, 4], [46, 4], [47, 4], [48, 4], [49, 4], [50, 4], [51, 4], [52, 4], [53, 4], [54, 4], [55, 4], [56, 4], [57, 5], [58, 5], [59, 5], [60, 5], [61, 4], [62, 4], [63, 4], [64, 4], [65, 5], [66, 5], [67, 5], [68, 5], [69, 5], [70, 5], [71, 5], [72, 4], [73, 5], [74, 5], [75, 5], [76, 5], [77, 5], [78, 5], [79, 5], [80, 5], [81, 5], [82, 6], [83, 6], [84, 6], [85, 5], [86, 5], [87, 5], [88, 5], [89, 5], [90, 5], [91, 6], [92, 6], [93, 6], [94, 6], [95, 6], [96, 6], [97, 6], [98, 5], [99, 5], [100, 5]]: plots:-display( [ plot(Points, style=point, color=red, symbol=solidcircle),                   plot(f, 10..100)                 ],                 gridlines=false               );
 > g:=n->ceil(sqrt(4*floor(n+1/2)))-floor(sqrt(2*floor(n+1/2)))-1:   plots:-display( [ plot(Points, style=point, color=red, symbol=solidcircle),                     plot(g, 10..100)                   ],                   gridlines=false                 );
 >

at

https://en.wikipedia.org/wiki/Rouch%C3%A9%E2%80%93Capelli_theorem

A system of linear equations with n variables has a solution if and only if the rank of its coefficient matrix A is equal to the rank of its augmented matrix [A|b].[1] If there are solutions, they form an affine subspace of ${\displaystyle \mathbb {R} ^{n}}$ of dimension n − rank(A). In particular:

1. if n = rank(A), the solution is unique,
2. otherwise there are infinitely many solutions.

It is trivial to produce the coefficient matrix, the augmented matrix and their ranks, using the code

```sys:=[subs(x = xmin, uhat1) = rhs(bcf1),
subs(x = xmax, uhat1) = rhs(bcf2),
subs(p = pmin, uhat1) = rhs(bcf3),
subs(p = pmax, uhat1) = rhs(bcf4)]:
vars:=[Af[0, 1](t), Af[0, 2](t), Af[0, 0](t), Af[0, 3](t)]:
coeffMat,b:=GenerateMatrix(sys, vars):
augMat:=GenerateMatrix(sys, vars, augmented=true):
Rank(coeffMat);
Rank(augMat);
```

which shows that the coefficient matrix has rank 3, and the augmented matrix has rank 4. Thus the sytem has no solution

Maple 2021

ScientificConstants[GetValue] is not working in at least Maple 2021 and Maple 2022

No mention of Maple Flow - if you had mentioned Maple Flow, I wouldn't even have attempted an answer

## Perhaps...

if you uploaded a usable worksheet someone might investigate this

## hmmm...

The first ode you see is a screen shot of the maple help page to show that the syntax using series worked. That was not an example of mine. It is just screen shot of the help page to show that the syntax in help page worked but gives error when I used it.

The first worksheet (odetest.mw) attached below is obtained directly from Maple help by using (in the help browser) the menu commands View -> Open Page as Worksheet. This whole worksheet has then been re-executed (using !!! in the Maple toolbar). Everything executes with no errors. So I still cannot replicate the problem implied by your remark

the syntax in help page worked but gives error when I used it.

And since you do not provide a worksheet illustrating the problem whihc you are experiencing - there isn't much I can do to fix it.

The second worksheet (testODE2.mw) below shows some of my attempts to replicate your problem. It generates series solutions for the same ODE from the help page in several slightly different ways, then tests each one with odetest() in a couple of different ways. Every single one of these "works" - so again I cannot replicate your problem

odetest

test the explicit or implicit results from ODE-solvers

 Calling Sequence Parameters Description Examples

 Calling Sequence odetest(sol, ODE, y(x)) odetest(sol, ODE, series, point=x0)

Parameters

 sol - Ordinary Differential Equation (ODE) solution being tested; can be a set or list of them ODE - ODE, or a set or list of them which can also include initial or boundary conditions y(x) - (optional) indeterminate function of the ODE or a set or list of them series - required when testing series solutions point=x0 - (optional) specification of the expansion point, x0, when testing series solutions

Description

 • The odetest command checks explicit and implicit solutions for ODEs by making a careful simplification of the ODE with respect to the given solution. If the solution is valid, the returned result will be 0; otherwise, the algebraic remaining expression will be returned. In the case of systems of ODEs, odetest can only test explicit solutions, given either as a set or as a list of sets. (For information on non-composed sets of solutions for nonlinear systems, see dsolve,system .)
 • To test whether a solution satisfies one or many initial or boundary conditions, pass to odetest the ODE together with the initial or boundary conditions, enclosed as a set  or list , as second argument.
 • If odetest returns a nonzero result, the solution being tested is not necessarily wrong; sometimes further simplifications or manipulations of odetest's output are required to obtain zero, and so verify the solution is correct. If the solution was obtained using the dsolve  command, it is recommended that you recompute the solution using one or both of the useInt and implicit options - see dsolve . This may facilitate the verification process. Also, an alternative testing technique, particularly useful with linear ODEs, is to try to recompute the ODE departing from the solution which odetest fails in testing. Examples of both types are found at the end of the next section.
 • To test series solutions , pass the keyword series as an extra argument. Only one series solution for one ODE (can be a set with initial/boundary conditions) can be tested.

Examples

An ODE problem with initial conditions

 > ODE := [diff(y(x),x)=sin(x-y(x)), y(0) = 8];
 (1)
 > sol := dsolve(ODE);
 (2)
 > odetest( sol, ODE );    # verifies 'sol' solves the ode and satisfies y(0) = 8
 (3)

A second order ODE problem with boundary conditions

 > ODE := diff(y(x),x,x) + diff(y(x),x) + y(x)=0;
 (4)
 > bc := y(0) = 1, y(2*Pi) = 1;
 (5)
 > sol := dsolve({ODE, bc});
 (6)
 > odetest( sol, [ODE,bc] ); # verifies 'sol' solves the ODE and satisfies the bc given
 (7)

A series solution for a nonlinear ODE with initial conditions

 > ODE := [diff(y(x),x,x)+diff(y(x),x)^2=0, y(a)=0, D(y)(a)=1];
 (8)
 > sol := dsolve( ODE, y(x), type='series');
 (9)
 > odetest(sol, ODE, series);
 (10)

When testing series solutions and the initial conditions are not present in the input to odetest, an indication of the expansion point is required

 > ODE := diff(diff(y(x),x),x) = (3*x^2+c)*diff(y(x),x)+((3-b)*x-a)*y(x);
 (11)
 > sol := y(x) = series(1+(-1/2*a)*x^2+(-1/6*b+1/2-1/6*c*a)*x^3+(1/24*a^2-1/24*c*b+1/8*c-1/24*c^2*a)*x^4+O(x^5),x,5);
 (12)
 > odetest(sol, ODE, series, point = 0);
 (13)

An ODE with an arbitrary function  of (x, y, dy/dx) and a solution involving nested integrals with a RootOf  in the integrand

 > ODE := diff(y(x),x,x) = 1/x^2*_F1(diff(y(x),x)*x/y(x))*y(x);
 (14)
 > sol := dsolve(ODE);
 (15)
 > odetest(sol,ODE);
 (16)

Testing ODE solutions given in implicit form, that is, not solved for the unknown (here y(x))

 > ODE := diff(y(x),x)=F((y(x)-x*ln(x))/x) + ln(x);
 (17)
 > sol := dsolve(ODE,implicit);
 (18)
 > odetest(sol,ODE);
 (19)

When the ODE has derivatives of other indeterminate functions and the solution is implicit, the specification of the indeterminate function of the problem is required by both dsolve  and odetest

 > ODE := diff(y(x),x) = x*f(x)^2*(x+2)*y(x)^3+f(x)*(x+3)*y(x)^2-diff(f(x),x)/f(x)*y(x);
 (20)
 > sol := dsolve(ODE, y(x), implicit);
 (21)
 > odetest(sol,ODE, y(x));
 (22)

Testing reductions of order returned by dsolve  using ODESolStructures

 > ODE := diff(y(x),x,x) = (diff(y(x),x)-y(x)^3-f(x)+3*x*y(x)^2*diff(y(x),x)+x*diff(f(x),x))/x;
 (23)
 > sol := dsolve(ODE, y(x));
 (24)
 > odetest(sol,ODE);
 (25)

A linear system of ODEs. The solution is a set containing x(t) and y(t) as functions of t.

 > sysODE := {diff(y(t),t)=-x(t),diff(x(t),t)=y(t)}, {x,y}(t);
 (26)
 > solsys := dsolve(sysODE);
 (27)
 > odetest(solsys,sysODE);
 (28)

A nonlinear system of ODEs. The solution is a list of sets, the first one containing the possible answers for x(t), and the second one expressing y(t) as a function of x(t):

 > sysODE := {diff(y(t),t)=-x(t)^2,diff(x(t),t)=y(t)}, {x,y}(t);
 (29)
 > solsys := dsolve(sysODE);
 (30)

These answers can be tests by passing them to odetest as a list.

 > odetest(solsys,sysODE);
 (31)

Alternatively, you can call dsolve  with the 'explicit' extra argument to directly obtain (many) composed solution sets. To test all these answers, use the map  function to apply odetest to each solution set:

 > solsys := dsolve(sysODE, explicit);
 (32)
 > map(odetest, [solsys], sysODE);
 (33)

One possible workaround for an example where odetest fails in verifying dsolve 's solution

 > ODE := diff(y(t),t) = ((b+2+2*t)*y(t)+1)/(1-(1+t)^2);
 (34)
 > sol := dsolve(ODE);
 (35)
 > odetest(sol, ODE);            # fails in verifying this solution
 (36)
 > sol2 := dsolve(ODE, useInt);  # compute 'sol' again with 'useInt'
 (37)
 > odetest(sol2, ODE);           # this solution is easier to test
 (38)

By evaluating the integrals appearing in sol2, the output returned by dsolve  without using the 'useInt' option can be constructed from the one obtained using the 'useInt' option, which was already verified to be correct.

An example hard to test due to the presence of radicals and Kummer  functions in the solution

 > ODE := diff(y(x),`\$`(x,2)) = ((-a[1]*F+a[0]*E)*x+B*e*a[1])*y(x)/B^2/e^2/x/E^3;
 (39)
 > sol := dsolve(ODE, [hyper3]);             # 1F1 and KummerU solution
 (40)
 > ode := PDEtools[dpolyform](sol, no_Fn);  # the ode satisfied by sol
 (41)
 > normal( ODE - op([1,1],ode) );             # ode = ODE
 (42)

Yet another alternative is to convert the special functions entering sol to other functions easier to test; in this example convert from Kummer  to Whittaker  functions:

 > sol_W:=convert(sol,Whittaker);
 (43)
 > odetest( sol_W, ODE );    # this Whittaker solution is easier to test
 (44)
 >

 > restart:   interface(version); # # ode from help page for odetest() command #   ODE := diff(diff(y(x),x),x) = (3*x^2+c)*diff(y(x),x)+((3-b)*x-a)*y(x);
 (1)
 > # # Solution from help page for odetest() command. # # NB this "solution" is given, not computed, and it # is 5-th order #   sol:= y(x) = series(1+(-1/2*a)*x^2+(-1/6*b+1/2-1/6*c*a)*x^3+(1/24*a^2-1/24*c*b+1/8*c-1/24*c^2*a)*x^4+O(x^5),x,5); # # Use dsolve to get a series solution for the above ODE # (no boundary conditions). This will be 6-th order (by # default) #   sol2_6:=  dsolve(ODE, y(x), series); # # Just for completeness (and comparison with sol1 above) # generate a 5-th order solution #   Order:=5:   sol2_5:=  dsolve(ODE, y(x), series); # # Use dsolve to get a series solution (with boundary conditions) # at both 5-th and 6-th order #   sol3_5:=  dsolve([ODE, y(0)=1, D(y)(0)=0], y(x), series);   Order:=6:   sol3_6:=  dsolve([ODE, y(0)=1, D(y)(0)=0], y(x), series)
 (2)
 > # # Five (slightly different) series solutions above. Check each # of them (two ways) #   odetest(sol, ODE, series, point = 0); # from maple help page   odetest(sol, ODE, 'series', point = 0);   odetest(sol2_5, ODE, series, point = 0);   odetest(sol2_5, ODE, 'series', point = 0);   odetest(sol2_6, ODE, series, point = 0);   odetest(sol2_6, ODE, 'series', point = 0);   odetest(sol3_5, [ODE, y(0)=1, D(y)(0)=0], 'series');   odetest(sol3_5, [ODE, y(0)=1, D(y)(0)=0], series);   odetest(sol3_6, [ODE, y(0)=1, D(y)(0)=0], 'series');   odetest(sol3_6, [ODE, y(0)=1, D(y)(0)=0], series);
 (3)
 >