## 4596 Reputation

9 years, 119 days

## Observation...

This is one of these case where (as Carl has said), it is probably easier to rewrite your ODE using brain-power rather than Maple. On the other hand, you can use the latter if you want to - see the attached

 > restart;   expr:= dx+4*x*dy/y = dy*a*y^(-6*b-1);   dsolve(algsubs(dy/dx=diff(y(x),x),subs(y=y(x),expand(expr/dx)) ));
 (1)
 >

## Like this?...

see the attached

 > restart;   c := x -> piecewise(0 <= x and x <= 450000, 0.37*x, 450000 < x, 0.37*x + 0.06*(x - 450000)):   g := x -> piecewise(0 <= x and x <= 558000, 0.37*x, 558000 < x, 0.37*x + 0.12*(x - 558000)):   plots:-animate( plot,                   [ [g(x), c(x)],                     x = 0..skat,                     color = [red, blue],                     legend=['g(x)', 'c(x)']                   ],                   frames = 20,                   skat = 0 .. 1000000                 );
 >
 >

## Errrr - works in Maple 2018!??...

see the attached

 > restart; interface(version); Physics:-Version(); sys:=[diff(Y(x, t), t, t) = 0, Y(x, 0) = 0, (D[2](Y))(x, 1) = 0]; pdsolve(sys);
 (1)
 >

## Use '~', the elementwise operator...

as in

```restart;
m1 := Matrix([[x[0,0]=1,x[0,1]=1,x[0,2]=0,x[0,3]=0,x[0,4]=0,x[0,5]=0],
[x[1,0]=0,x[1,1]=0,x[1,2]=1,x[1,3]=1,x[1,4]=0,x[1,5]=0],
[x[2,0]=0,x[2,1]=0,x[2,2]=0,x[2,3]=0,x[2,4]=1,x[2,5]=1]
]
);
m2:=rhs~(m1);
```

## Requirement no completely clear...

but the attached *may* be what you want

 > expr:=u[1](T[1], T[2], T[3])=exp(T[1]*I); diff(expr, indets(rhs(expr),name)[]);
 (1)
 >

## One way...

The method shown in the attached is pretty 'ugly', but it seems to work

 > restart;   getOrd:= (f1, f2)-> max                       ( PDETools:-difforder                                   ( select                                     ( i-> has(i, diff(f2,t)),                                       convert                                       ( indets(f1, function),                                         list                                       )                                     )                                   )                       ):
 > f:= 2*y*(diff(x(t), t))^2+3*(diff(x(t), t\$3))-3*x*(diff(y(t), t));   getOrd(f, y(t));   getOrd(f, x(t));
 (1)
 >

## Probably several ways to do this...

some of which are probably more "elegant" than my attempt in the attached

 > restart;   plot3d( sin(x*y),           x = -2 .. 2,           y = -2 .. 2,           axes = framed,           colorscheme = [ "zgradient",                           ["Blue", "Green", "Yellow", "Red"]                         ]         );
 > doContour:= proc( nsteps::posint)                     local j;                     uses plots, ColorTools:                     return display                            ( [ seq                                ( implicitplot                                  ( sin(x*y)=-1+2*(j-1)/nsteps,                                    x = -2 .. 2,                                    y = -2 .. 2,                                    axes = framed,                                    color = WavelengthToColor                                            ( 450+(j-1)*(625-450)/nsteps,                                              method="linear"                                            ),                                    if   type(j, odd)                                    then legend=(-1+2*(j-1)/nsteps)                                    else NULL                                    fi,                                    scaling=constrained,                                    size=[800,800]                                  ),                                  j = 1..nsteps+1                                )                              ]                            );               end proc:  doContour(10);  doContour(20);  doContour(100);
 >
 >
 >

## Format...

of your animate() command is a bit "unusual". In fact I was mildly urprsed that it "worked" at all!

If I convert it to follow the standard technique explained in the "help", then the line properties appear to be retained. See the attached

 > with(plots): A := animate( plot,               [ X+t,                 X = -10 .. 10,                 color = red,                 linestyle = dash               ],               t = -10 .. 10             ): B := animate( plot,               [ t^2+X,                 X = -10 .. 10,                 color = blue,                 linestyle = solid               ],               t = -10 .. 10             ): display([A, B]);
 >

## So many things to fix...

that it is impossible to list them all.

And so many that it is also impossible (for me at least) to figure out what you actually want.

See the attached for details of basix syntax (and other) problems. Be aware I got bored and stopped identifying errors about half way through your worksheet

 >
 >
 >
 >
 >
 >
 >

## A fundamental problem...

Outputting solution values to a file (or multiple files) is trivial.

However I'm not going to show how to do this because I have conerns over the validity of any generated numerical values.

Consider very carefully your boundary conditions

cond:= {u(0,t)=sin(.5*Pi), u(x,0)=0}

The first condition states that                u(0,0)=sin(0.5*Pi)=1.
The second condition states that          u(0,0)=0

So what do you want u(0,0) to be? Produce consistent  boundary conditions accordingly

## Hmmm...

For any given value of 's', you expresision is quartic in 'z' - so has four roots

A quick check examination appears to show that for any(?) value of 's', the expression has two real roots and a conjugate pair for 'z'.

The attached shows the imaginary parts of the conjugate pair

 > restart;   expr:=.9215999996*z^4+.9999999996*z^4*s^4-2.079999999*z^4*s^2-3.212799999*s*z^3+12.63999999*s^3*z^3-7.999999997*s^5*z^3-24.17694443*s^4*z^2-.3071555554*s^2*z^2+23.99999999*s^6*z^2-31.99999999*s^7*z+6.804911108*s^3*z+11.58777777*s^5*z+15.99999999*s^8+5.692222220*s^6-4.853219442*s^4=0:   f:=p-> min(Im~([fsolve(eval(expr,s=p), z, complex)])):   g:=p-> max(Im~([fsolve(eval(expr,s=p), z, complex)])):   plot( [f, g], -10..10, axes=boxed);
 >

## The procedure dicho()...

has two problems (at least)

I have highlighted both of these in the attached

 > restart;   with(DEtools):   eq1:= (D(x))(t) = -y(t);   eq2:= (D(y))(t) = x(t)+2*x(t)^3-signum(z(t));   eq3:= (D(z))(t) = w(t);   eq4:= (D(w))(t) = -z(t)*(1+6*x(t)^2);   sys:= eq1, eq2, eq3, eq4;   ic1:= [ x(0)=0, y(0)=0, z(0)=cos(1), w(0)=sin(1)];   ic2:= [ x(0)=0, y(0)=0, z(0)=cos(2.5), w(0)=sin(2.5)];   ic:= ic1, ic2;   DEplot( [sys],           [x(t), y(t), z(t), w(t)],            t = 0 .. 10,            [ic],            stepsize = 0.5e-1,            scene = [x(t), y(t)],            linecolor = [blue, red]         );   sol1:= dsolve          ( { sys, x(0)=0, y(0)=0, z(0)=cos(1), w(0)=sin(1)},            { x(t), y(t), z(t), w(t)},            type=numeric          );   T:= 10.0;   N:= 100;   h:= T/N;   xk:= 0;   for k from 1 to N do       solk:= sol1(k*h);       xknew:= subs(solk,x(t));       yknew:= subs(solk,y(t));       if   xk*xknew<=0 and abs(yknew-6)<0.5       then break       fi;   xk:= xknew;   od:   sol1(k*h):   temps:= proc(alpha,eps)                 local sol, solk, T, N, h ,k, xk, xknew,                       yknew, t0, t1, tm, x0, x1, xm;                 sol:= dsolve                        ( { sys, x(0)=0, y(0)=0, z(0)=cos(alpha) ,w(0)=sin(alpha) },                          { x(t), y(t), z(t), w(t)},                          type=numeric                        );                 T:= 10.0;                 N:= 100;                 h:= T/N;                 xk:= 0;                 for k from 1 to N do                     solk:= sol(k*h);                     xknew:= subs(solk, x(t));                     yknew:= subs(solk, y(t));                     if   xk*xknew<=0 and abs(yknew-6)<0.5                     then break                     fi;                     xk:= xknew;                 od;                 t0:= (k-1)*h;                 t1:= k*h;                 x0:= subs(sol(t0), x(t));                 x1:= subs(sol(t1), x(t));                 while abs(x0-x1)>eps do                       tm:= (t0+t1)/2;                       xm:= subs(sol(tm), x(t));                       if   xm*x0<0                       then x1:= xm;                            t1:= tm;                       else x0:= xm;                            t0:= tm;                       fi;                 od;                 return t0;            end:   dicho:= proc(eps)                 local a, b, m, sola, solb, solm,                       ta, tb, tm, ya, yb, ym;                 a:= 1;                 b:= 2.5;                 sola:= dsolve                         ( {sys, x(0)=0, y(0)=0, z(0)=cos(a), w(0)=sin(a)},                           {x(t),y(t),z(t),w(t)},                           type=numeric                         );                 solb:= dsolve                         ( {sys, x(0)=0, y(0)=0, z(0)=cos(b), w(0)=sin(b)},                           {x(t), y(t), z(t), w(t)},                           type=numeric                         );                 ta:= temps(a,eps);                 tb:= temps(b,eps);                 ya:= subs(sola(ta),y(t));                 yb:= subs(solb(tb),y(t));                 while abs(yb-ya)>eps do                       m:= evalf((a+b)/2);                       solm:= dsolve                               ( { sys, x(0)=0, y(0)=0, z(0)=cos(m), w(0)=sin(m)},                                 {x(t), y(t), z(t), w(t)},                                 type=numeric                               );                       tm:= temps(m,eps);                       # yb:= subs(sol(tm),y(t)); # should this be sola, solb or                                                  # solm?? Used solm() just to                                                  # eliminate the error This is                                                  # quite likely to be incorrect                       yb:= subs(solm(tm),y(t));                       if   (ym-6)*(ya-6)<0 # first time through this loop ym is                                            # not initialized In fact it is                                            # *never* initialized - although all                                            # sorts of stuff depends on it.                                            # WTF is ym - I have no idea!!                       then b:= m;                            yb:= ym;                       else a:= m;                            ya:= ym;                       fi;                 od;                 return a;             end:   dicho(0.01);   temps(2.136718750,0.01);
 (1)
 >

## If...

I understand correctly, then you want something like the attached

 > restart; # # Define the functions # #      u(t,f) and #      x(t,f) #   u:= unapply       ( 1-(8*(10.3968*t^2-5.8368*t*f-.229376*f^2-5.1984))/(4.56*t^2-2.56*t*f+.8192*f^2+2.28)^2,         [t,f]       ):   x:= unapply       ( f+t-(8*(-2.28*t+.64*f))/(2.28+2*(-t+.64*f)^2+2.56*t^2),         [t,f]        ): # # Plot u(0,f) against x(0,f) #   plot( [x(0,f), u(0,f), f=-5..5]); # # Plot u(t,f) against x(t,f) and t #   plot3d([ x(t,f), t, u(t,f) ], f=-5..5, t=-5..5);
 >

## Other this reasons this procedure will n...

For details, see the attached.

Short version is that your procedure calls two functions 'f()' and 'check()' which are defined nowhere so procedure is guaranteed to fail, unless these 'unidentified' functions exist elsewhere(probably top-level)  in your code.

The "debug" icon in the toolbar is only really intended for to work with currenlty running operations - so if some code seems to be taking "forever" to execute then this will give you a kind of "interrupt" whihmay help you to figure out what it is currently doing and why. I actually find this facility pretty much useless for code which fails "quickly".

The best way to work out how to use the debugger (beware it isn't very user-friendly!) is to read the Maple Programming Guide available here, and experiment with the examples

https://www.maplesoft.com/support/help/Maple/view.aspx?path=ProgrammingGuide/Chapter16

(This document used to be available as a pdf (easier to download and study, in my opinion) but I can' figure out where the pdf link is

 > f1:=  proc(n)           local x, y, i, rx, ry, dumbindex:= 1,                 xdumb, ydumb, d:= 0.00197241, flag:= 1,                 p, xl, yl, temp, xg, yg;           xdumb[1]:= n:           ydumb[1]:= n:             for i from 1 to  n do                   while (d>0.0019724 or flag=1) do                   #                   # Why is a new random number generator                   # being instantiated every time this                   # loop executes?                   # Surely better would be to initialise                   # the first RNG outside both loops, call                   # and run the calls to it in the outer                   # - or maybe not?? Even if a couple of                   # 'new' random numbers are needed every                   # the 'inner' loop runs, it is probably                   # a better idea to inilis the generators                   # in the outer loop and only have the                   # 'calls' here                   #                        rx:= rand(1..n):                     x[i]:= rx():                     ry:= rand(x[i]..n):                     y[i]:= ry():                   #                   # What does the following command do??                   # There is a command called 'check()'                   # in the Physics() package, but since                   # the latter has not been loaded(?)                   # the check() call here will never return                   # anything other than an unevaluated                   # function call??                   #                   # Doesn't look like right!                   #                     flag:= check(x[i],y[i],xdumb,ydumb,dumbindex):                       if   flag=1                     then dumbindex:= dumbindex+1;                          xdumb[dumbindex]:= x[i]:                          ydumb[dumbindex]:= y[i]:                     fi:                   #                   # What does the next command do??                   #                   # The function 'f()' is undefined so will                   # never return anything other than its                   # literal name, which is a single value.                   # Therefore no way that a duplicate                   # assignment will work, because lhs has                   # two names, 'p' and 'd', and rhs only has                   # one (unevaluated) value.                   #                   # Since both 'p' and 'd' are used in this                   # code - this statement is probably not a                   # good idea                   #                   # I did consider that this *might* be a                   # recursive call - ie calling this procedure                   # with new arguments, but since this                   # procedure is called 'f1' and accepts one                   # argument, and the invocation calls 'f'                   # and supplies two arguments this interpretation                   # is almost certainly incorrect                   #                   # The mismatch between the returned                   # (unevaluated) single value and the dual                   # assignment is what causes the error message                   # whenever the prodedure 'f1()' is invoked                   #                      p, d:= f(x[i],y[i]):               end:                 xl[i]:= x[i]:               yl[i]:= y[i]:                 if   i=1               then temp:= p:                    xg:= x[i]:                    yg:= y[i];                 else if   p
 (1)
 > f1(10);
 >

## There are a few typos...

whihc are pretty simple to fix. However a more serious problem is in the form of the boundary conditions which you are using. You are only allowed derivatives which are "normal" to the boundary. You have two,

D[2](C)(0, t) = 0, D[2](C)(250, t) = 0

which are not! These can never be used

An alternative is that you have misinterpreted how to apply the operator 'D'. The two boundary conditions above mean

1. Obtain the first derivative with respct to the second variable (ie 't')
2. Evaluate the results of (1) at (0,t) and (250,t)

On the other hand if you had intended to obtain the second derivative with respect to the first variable, this would be written

D[1,1](C)(0, t) = 0, D[1,1](C)(250, t) = 0

These derivatives are "normal" to the boundary and so pdsolve/numeric will provide a solution

I have used this changed interpretation in the final two execution groups of the attached, which now works without error. However a quick check on a few values indicates that the result isn't numerically very stable. When t>~0.15 values start to grow very  rapidly

 >
 >
 >
 >
 >
 >
 >
 >
 >
 (1)
 >
 >
 >

 1 2 3 4 5 6 7 Last Page 1 of 105
﻿