tomleslie

12497 Reputation

19 Badges

12 years, 231 days

MaplePrimes Activity


These are replies submitted by tomleslie

@jud 

I get the same error as you do in Maple 2021. Looks like a bug that has been fixed

@ijuptilk 

why do I have repeated values for j and u? And how can solve this problem without using Makeunique function.

Thhis is because the loop system I set up treats 'u' and 'j' as independent variables, so in the worksheet I provided initially, the command

 ans1:= CodeTools:-Usage(Array([seq(seq( [j, doCalc(j,u)], j=-2..0, 2), u =1..334,333)])):

sets u=1, then does the calculation for j=-2 and 0, the sets u=334 amd does the calculation for j=-2 and 0. It is still not clear to me what you actually want to achieve, unless 'j' is not an independent variable, but is in fact a simple function 'u'. One possibility is shown in the attached where as u varies from 1..334, j varies from -2..0, according to j=-2.0+(u-1)/334.

(In the attached, purely to save execution time I do this calculation with 'u' varying in steps of 33. If you want to do this caclulation ffor all u-values from 1..334, then just change 'ustep' to 1

#######################################################

I also want to do contourplot for u vs j for complex values of p_min. 

So far as I can tell all returned "values" of p_min are expressions containing the variables 'x' and 't'. Maybe some of these expression will be complex for somevlaues of 'x' and 't', and maybe not. Since I have absoutely no idea what values of 'x' and 't' you paln on using, I haven't really considered this question further


 

  restart:
  doCalc:= proc( xi , u)  #u is the \bar(H): normalize magnetic field magnitude,
                          # where H = bar(H)*H__a
                 # Import Packages
                 uses ArrayTools, Student:-Calculus1, LinearAlgebra,
                      ListTools, RootFinding, plots, ListTools:
                 local gamma__1:= .1093,
                       alpha__3:= -0.1104e-2,
                       k__1:= 6*10^(-12),
                       d:= 0.2e-3,
                       theta0:= 0.1e-3,
                       eta__1:= 0.240e-1, chi:= 1.219*10^(-6),
                       alpha:= 1-alpha__3^2/(gamma__1*eta__1),
                       theta_init:= theta0*sin(Pi*z/d),
                       H__a:= Pi*sqrt(k__1/chi)/d,
                       H := u*H__a,     
                       c:=alpha__3*xi*alpha/(eta__1*(4*k__1*q^2/d^2-alpha__3*xi/eta__1 - chi*H^2)),
                       w := chi*H^2*eta__1*alpha/(4*k__1*q^2/d^2-alpha__3*xi/eta__1 - chi*H^2),
                       n:= 20,
                       g, f, b1, b2, qstar, OddAsymptotes, ModifiedOddAsym,
                       qstarTemporary, indexOfqstar2, qstar2, AreThereComplexRoots,
                       soln1, soln2, qcomplex1, qcomplex2, gg, qq, m, pp, j, i,
                       AllAsymptotes, p, Efun, b, aa, F, A, B, Ainv, r, theta_sol, v, Vfun, v_sol,minp,nstar;

# Assign g for q and plot g, Set q as a complex and Compute the Special Asymptotes

                 g:= q-(1-alpha)*tan(q)- (w*q + c)*tan(q):
                 f:= subs(q = x+I*y, g):
                 b1:= evalc(Re(f)) = 0:
                 b2:= evalc(Im(f)) = 0:
                 qstar:= (fsolve(1/c = 0, q = 0 .. infinity)):
                 OddAsymptotes:= Vector[row]([seq(evalf(1/2*(2*j + 1)*Pi), j = 0 .. n)]);

# Compute Odd asymptote

                 ModifiedOddAsym:= abs(`-`~(OddAsymptotes, qstar));
                 qstarTemporary:= min(ModifiedOddAsym);
                 indexOfqstar2:= SearchAll(qstarTemporary, ModifiedOddAsym);
                 qstar2:= OddAsymptotes(indexOfqstar2);

# Compute Odd asymptote

                 AreThereComplexRoots:= type(true, 'truefalse');
                 try
                      soln1:= fsolve({b1, b2}, {x = min(qstar2, qstar) .. max(qstar2, qstar), y = 0 .. infinity});
                      soln2:= fsolve({b1, b2}, {x = min(qstar2, qstar) .. max(qstar2, qstar), y = -infinity .. 0});
                      qcomplex1:= subs(soln1, x+I*y);
                      qcomplex2:= subs(soln2, x+I*y);
                 catch:
                       AreThereComplexRoots:= type(FAIL, 'truefalse');
                 end try;

# Compute the rest of the Roots

                 OddAsymptotes:= Vector[row]([seq(evalf((1/2)*(2*j+1)*Pi), j = 0 .. n)]);
                 AllAsymptotes:= sort(Vector[row]([OddAsymptotes, qstar]));
                 if   AreThereComplexRoots
                 then gg:= [ qcomplex1,
                             qcomplex2,
                             op(Roots(g, q = 0.1e-3 .. AllAsymptotes[1])),
                             seq(op(Roots(g, q = AllAsymptotes[i] .. AllAsymptotes[i+1])), i = 1 .. n)
                          ];
                 elif not AreThereComplexRoots
                 then gg:= [ op(Roots(g, q = 0.1e-3 .. AllAsymptotes[1])),
                             seq(op(Roots(g, q = AllAsymptotes[i] .. AllAsymptotes[i+1])), i = 1 .. n)
                           ];
                 end if:

# Remove the repeated roots if any & Redefine n

                 qq:= MakeUnique(gg):
                 m:= numelems(qq):

## Compute the `τ_n`time constants

                 for i to m do
                 p[i] := gamma__1*alpha/(4*k__1*qq[i]^2/d^2 - alpha__3*xi/eta__1- chi*H^2):
                 end do:

## Compute θ_n from initial conditions

                for i to m do
                Efun[i] := cos(qq[i]-2*qq[i]*z/d)-cos(qq[i]):
                end do:

## Compute integral coefficients for off-diagonal elements θ_n matrix

                printlevel := 2:
                for i to m do
                    for j from i+1 to m do
                        b[i, j] := int(Efun[i]*Efun[j], z = 0 .. d):
                        b[j, i] := b[i, j]:
                        aa[i, j] := b[i, j]:
                        aa[j, i] := b[j, i]:
                    end do :
                end do:

## Calculate integral coefficients for diagonal elements in theta_n matrix

                for i to m do
                aa[i, i] := int(Efun[i]^2, z = 0 .. d):
                end do:

## Calculate integrals for RHS vector

               for i to m do
               F[i] := int(theta_init*Efun[i], z = 0 .. d):
               end do:

## Create matrix A and RHS vector B

               A := Matrix([seq([seq(aa[i, j], i = 1 .. m)], j = 1 .. m)]):
               B := Vector([seq(F[i], i = 1 .. m)]):

## Calculate inverse of A and solve A*x=B

              Ainv := 1/A:
              r := MatrixVectorMultiply(Ainv, B):

## Define Theta(z,t)
             theta_sol := add(r[i]*Efun[i]*exp(-t/p[i]), i = 1 .. m):

## Compute v_n for n times constant

             for i to m do
             v[i] := (-2*k__1*alpha__3*qq[i])*(1/(d^2*eta__1*alpha*gamma__1))+ alpha__3^2*xi/(2*(eta__1)^2*qq[i]*alpha*gamma__1)+xi/(2*eta__1*qq[i]) + alpha__3*chi*H^2/(2*eta__1*qq[i]*gamma__1*alpha):
             end do:

## Compute v(z,t) from initial conditions
             for i to m do
             Vfun[i] := d*sin(qq[i]-2*qq[i]*z/d)+(2*z-d)*sin(qq[i]):
             end do:

## Define v(z,t)
             v_sol := add(v[i]*r[i]*Vfun[i]*exp(-t/p[i]), i = 1 .. m):

##
             minp:=min( seq( Re(p[i]), i=1..m) ):
             member(min(seq( Re(p[i]), i=1..m)),[seq( Re(p[i]), i=1..m)],'nstar'):

## Return all the plots
                 return theta_init, theta_sol, v_sol, minp, eval(p), nstar, p[nstar], g, H, H__a;
                 end proc:

# Run the calculation for supplied value of 'xi'

# Put the returned quantities in a simple list

ans:=[doCalc(-0.06, 0.001)]:
ans1:=[doCalc(0.06, 0.001)]:
evalf(ans[9]);
evalf(ans[10]);

0.3484926715e-1

 

34.84926715

(1)

##

with(plots):
d:= 0.2e-3:

## director Plot negative activity
display( plot(ans[1], z=0..d, color = green),
         plot([seq(eval(ans[2], t = j), j in [seq(.1*j, j = 1 .. 5)])], z=0..d, legendstyle = [font                 = ["ROMAN", 12],location = right], labeldirections = ["horizontal", "vertical"]),                     axis[1]=[tickmarks=[0="0", 0.00005="50", 0.00010="100", 0.00015="150", 0.00020=                         "200"]], axis[2]=[tickmarks=[0="0", 0.00002="2e-05", 0.00004="4e-05", 0.00006=                     "6e-05", 0.00008= "8e-05", 0.0001="1e-04"]], axes=boxed, labels =[typeset(z," (`μm`)                     "), typeset(theta(z, t)," (degrees)")], labelfont = ["TimesNewRoman", 14],                             axesfont = [14, 14]);


## director angle Plot postive activity
display( plot(ans1[1], z=0..d, color = green, legend = theta[0]),
         plot([seq(eval(ans1[2], t = j), j in [seq(.1*j, j = 1 .. 5)])], z=0..d, legend = [t[1] =                  .2, t[2] = .4, t[3] = .6, t[4] = .8, t[5] = 1.0], legendstyle = [font = ["ROMAN", 12],                location = right], labeldirections = ["horizontal", "vertical"]), axis[1]=[tickmarks=                [0="0", 0.00005="50", 0.00010="100", 0.00015="150", 0.00020="200"]], axis[2]=                        [tickmarks=[0="0", 0.00002="2e-05", 0.00004="4e-05", 0.00006="6e-05", 0.00008=                       "8e-05", 0.0001="1e-04"]], axes=boxed, labels =[typeset(z," (`μm`)"), typeset(theta(z, t)," (degrees)")], labelfont =["TimesNewRoman", 14], axesfont = [14, 14]);

## v Plot negative activity
plt2:= plot([seq(eval(ans[3], t = j), j in [seq(.2*j, j = 1 .. 5)])], z = 0 .. d, legendstyle = [font = ["ROMAN", 12],                location = right], labeldirections = ["horizontal", "vertical"]):

display( plt2, axis[1]=[tickmarks=[0="0", 0.00005="50", 0.00010="100", 0.00015="150", 0.00020="200"]          ], axes=boxed, labels =[typeset(z," (`μm`)"), typeset(v(z, t)," (m/s)")], labelfont = ["TimesNewRoman", 14], axesfont = [14, 14]);
 
## v Plot postive activity
plt21:= plot([seq(eval(ans1[3], t = j), j in [seq(.2*j, j = 1 .. 5)])], z = 0 .. d, legend = [t[1] =                .2, t[2] = .4, t[3] = .6, t[4] = .8, t[5] = 1.0], legendstyle = [font = ["ROMAN", 12],                location = right], labeldirections = ["horizontal", "vertical"]):

display( plt21,axis[1]=[tickmarks=[0="0", 0.00005="50", 0.00010="100", 0.00015="150", 0.00020="200"]          ], axes=boxed, labels =[typeset(z," (`μm`)"), typeset(v(z, t)," (m/s)")], labelfont = ["TimesNewRoman", 14], axesfont = [14, 14]);

  

Error, missing operator or `;`

 

# q plot

#

plot(ans[8], q=0..5*Pi,view=[DEFAULT,-10..10], labels =[q, r(q)], labelfont = ["TimesNewRoman", 14], labeldirections = ["horizontal", "vertical"] );

 

# p values

evalf(ans[5]);

table( [( 1 ) = -13.88378137-65.82711280*I, ( 2 ) = -13.88378137+65.82711280*I, ( 3 ) = 10.67051633, ( 4 ) = 3.198773762, ( 5 ) = 1.566666433, ( 6 ) = .9327425647, ( 7 ) = .6194882059, ( 9 ) = .3307289368, ( 8 ) = .4415527478, ( 11 ) = .2054822445, ( 10 ) = .2570092031, ( 13 ) = .1399934044, ( 12 ) = .1680476193, ( 15 ) = .1014876810, ( 14 ) = .1184258882, ( 18 ) = 0.6788038186e-1, ( 19 ) = 0.6033276905e-1, ( 16 ) = 0.8794198385e-1, ( 17 ) = 0.7693928692e-1, ( 20 ) = 0.5397792362e-1, ( 21 ) = 0.4857705137e-1 ] )

(2)

 # Director angle plot for at d/2 as a funtion of time

theta_f := subs(z=d/2, Re(ans[2])):
plot(theta_f, t=0..1000, view=[DEFAULT, DEFAULT], labels =[typeset(t," (seconds)"), typeset(theta('d/2', t),"(rad)")], labelfont = ["TimesNewRoman", 14], axesfont = [14, 14], labeldirections = ["horizontal", "vertical"], axes=boxed);

theta_f := subs(z=d/2, ans1[2]):
plot(theta_f, t=0..1, NULLview=[DEFAULT,DEFAULT], labels =[typeset(t," (seconds)"), typeset(theta('d/2', t),"(rad)")], labelfont = ["TimesNewRoman", 14], axesfont = [14, 14], labeldirections = ["horizontal", "vertical"], axes=boxed);
 

 

 

#Convert and store xi and doCalc into an array
  #ans1:= Array([seq( [j, doCalc(j)], j=-2..0, 0.006)]):
 ustep:=33:
 ans1:= CodeTools:-Usage(Array([seq([-2.0+(u-1)/334, u, doCalc(-2.0+(u-1)/334, u)], u =1..334, ustep)])):

memory used=4.32GiB, alloc change=32.00MiB, cpu time=47.36s, real time=43.19s, gc time=5.96s

 

#

#

p_min:= convert(ans1[..,5],list):
indets~(p_min, name);
j_value := convert(ans1[..,1],list);
u_value := convert(ans1[..,2],list);
# plot( convert(ans1[..,1],list), p_min/30)
eval(p_min[1], [z=1, t=1]);

[{t, z}, {t, z}, {t, z}, {t, z}, {t, z}, {t, z}, {t, z}, {t, z}, {t, z}, {t, z}, {t, z}]

 

[-2.000000000, -1.901197605, -1.802395210, -1.703592814, -1.604790419, -1.505988024, -1.407185629, -1.308383234, -1.209580838, -1.110778443, -1.011976048]

 

[1, 34, 67, 100, 133, 166, 199, 232, 265, 298, 331]

 

HFloat(-0.006123129246163778)

(3)

## please why do I have repeated values for j and u? And how can solve this problem without using Makeunique function.

## I also want to do contourplot for u vs j for complex values of p_min.

NULL

Download sProc2.mw

You further state

 

This appears to work for me, see the attached. It may be a version issue: I'm running Maple 2022, wh8hc Maple version are you running?

  restart:
  interface(version);
  with(GroupTheory):
  CompositionSeries(SmallGroup(336,209));

`Standard Worksheet Interface, Maple 2022.0, Windows 7, March 8 2022 Build ID 1599809`

 

_m880103392

(1)

 

Download GT.mw

because obviously the attached, using Maple's inbult numeric solvers, just *works"

So why do anything else?

However if you do (for some strange reason)  wabt to use the "differential transformation method", then I would mak the following observations

  1. I know nothing about the differential transformation method
  2. It appears to be related to "power series" solution methods. In fact there appears to be some controversy, in the online literature  as to whether it is actually different from a Taylor series method, or is simply a disguised variant.
  3. The one thing that *all power series*  methods share is that they rely on a power series expansion about a single point. So they are appropriate for initial-value problems, where all  (so-called) boundary conditions refer to a single point - in other word they are not "boundary conditions", they are in fact "initial conditions"
  4. Your problem has boundary conditions defined at eta=0 and eta=infinity - so about which of these points do you wish to perform a power series expansion?
  5. I have no doubt that some here will suggest that I convert the OP's  boundary-value problem to an intial-value problem using the shooting method. Definitely possible, but I am seriously reluctant to do this when one of the boundaries is "infinity"

  restart:
  eqn1 := diff(f(eta), `$`(eta, 3))+f(eta)*(diff(f(eta), `$`(eta, 2)))-(diff(f(eta), eta))^2-lambda*(diff(f(eta), eta)) = 0:
  eqn2 := diff(theta(eta), `$`(eta, 2))+f(eta)*(diff(theta(eta), eta))*Pr = 0:
  Bcs := (D(f))(0) = 1, f(0) = 0, (D(f))(infinity) = 0, theta(0) = 1, theta(infinity) = 0:
  params:=[lambda = .5, Pr = 6.3, infinity=10]:
  sol:=dsolve( eval([eqn1, eqn2, Bcs], params), numeric):
  plots:-odeplot( sol,
                  [ [ eta, f(eta) ],
                    [ eta, theta(eta) ]
                  ],
                  eta=0..10,
                  color=[red, blue]
                );

 

 


 

Download odeSol.mw

@imparter 

misunderstanding of viable boundary conditions. Thje code

with(DETools):
with(plots):
with(IntegrationTools):
Nb:=1:Nt:=1:h(x):=exp(x):
Eq1 := (diff(f(x),x,x))+(((1/x)*(diff(f(x),x))+Nb*((diff(f(x),x))*diff(g(x),x))+Nt(diff(f(x),x)^2))):   
Eq2 := (diff(g(x),x,x))+(((1/x)*(diff(g(x),x))+(Nb/Nt)*((diff(f(x),x,x))+(1/x)*diff(f(x),x)))):       

# Boundary condition f(x)=0, g(x)=0 at x=h(x)" where h(x) is some function of x", 

 #f '(x)=0, g'(x)=0 at x=0

Cd1 := f(h(x)) = 0, (D(f))(0) = 0 :
dsys := {Cd1, Eq1}:
dsol := dsolve(dsys, numeric, output = operator):
plots[odeplot](dsol, [x, diff(f(x), x$1)], 0 .. 5, color = green):
Cd2 := g(h(x)) = 0, (D(g))(0) = 0:
dsys := {Cd1, Cd2, Eq1, Eq2}:
dsol := dsolve(dsys, numeric, output = operator):
plots[odeplot](dsol, [x, f(x)], 0 .. 5, color = red);
plots[odeplot](dsol, [x,g(x)], 0 .. 5, color = black);

contains two *interesting* definitions, namely

h(x):=exp(x)
f(h(x))=0

So just consider these for a moment. For any (real) value of x from -infinity to +infinity, exp(x) varies from 0 to +infinity, so h(x) is any value from 0 to +infinity - and the you define f(h(x))=0 - so f() for any value from 0 to +infinity is zero. You think this is a valid "boundary" condition??

@falkner 

there is no "tutor" available for use with dsolve().

One of the stategies used by  the dsolve() command is to try to identify whether a supplied ODE is of a certain "type".

The code

#
# Define system and solve
#
  sys:=[ diff(y(x),x) = y(x)^4 * cos(x) + y(x)*tan(x),
         y(0) = 0.5
       ]:
  DETools:-odeadvisor(sys[1]);
                          [_Bernoulli]

shows that Maple can identify the ODE as a Bernoulli equation. If you change the odeadviser() command in the above to

DETools:-odeadvisor(sys[1], help);

then Maple will display  help page describing the general strategy used for solving equations of this type

@Anthrazit 

there is very often more than one way to do something

is to always upload a worksheet which illustrates your problem!

Use the big green up-arrow in the Mapleprimes toolbar

@imparter 

to address the basic problem.

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));

{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)}

(1)

  DE2 := series(eval(de2, g = G), p = 0, n + 1):
  indets(DE2, function(name));

{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)}

(2)

``

Download indets.mw

@imparter 

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)

@imparter 

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
#

Error, (in dsolve) found the following equations not depending on the unknowns of the input system: {1}

 

#
# 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)});

b[0](x) = (1/10)*x^2

(1)

 

Download diffODEs.mw

 

 

@FDS 

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);

Vector(9, {(1) = 124.0, (2) = 130.0, (3) = 130.0, (4) = 119.0, (5) = 136.0, (6) = 118.0, (7) = undefined, (8) = 130.0, (9) = 95.0})

 

HFloat(122.75)

(1)

 

Download getMean.mw

@Muhammad Usman 

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
_______________________________________________________________________________

 

 

Download printStuff.mw

@ 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:-(

@JAMET 

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);

[2000, [10]], [2001, [4, 7]], [2002, [9, 12]], [2003, [6]], [2004, [2, 8]], [2005, [5]], [2006, [1, 10]], [2007, [4, 7]], [2008, [6]], [2009, [2, 3, 11]], [2010, [8]], [2011, [5]], [2012, [1, 4, 7]], [2013, [9, 12]], [2014, [6]], [2015, [2, 3, 11]], [2016, [5]], [2017, [1, 10]], [2018, [4, 7]], [2019, [9, 12]], [2020, [3, 11]], [2021, [8]], [2022, [5]]

(1)

 

 

Download fri_2.mw

 

 

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