tomleslie

5194 Reputation

15 Badges

9 years, 244 days

MaplePrimes Activity


These are replies submitted by tomleslie

which I probably shoould have made in my earlier answer.

If the grid points were cacluated using exact arithmetic (ie rational fractions), and the field strength computed using these "exact" values, before being converted to floating-point form, this problem should not(?) arise.

@Adam Ledger 

For whatever reason, when you issue the ListDirectory() command, it returns a list of filenames. It appears that some of these contain non-printing, non-ASCII characters. In order to manipulate these filenames in Maple, you are first going to have to get rid of these "weird" characters. The command

Implode(select( IsASCII, Explode(a)));;

will do this. If the input string contains only ASCII characters, then the command does nothing at all, just returns the original string. I've checked my original worksheet in MAple 18, Maple 2015, Maple 2016, Maple 2017, Maple 2018 and Maple 2019 - it works in all of them. So unless yu can demonstrate otherwise, there is no need to worry about the "build" version

 

@Carl Love 

Confused "strongly connected" with "complete". Must remember to check my graph theory terminology

The betweenness centrality of a vertex (say 'v') is essentially a measure of how often 'v' appears in the shortest paths between two other vertexes in the graph. However, in an unweighted strongly connected graph, the shortest path between any two vertexes is always 1 (with no "intermediate" vertexes).

Hence the betweenness centrality for every vertex in an unweighted strongly connected graph is zero

Or am I missing something??

@LasseFrost 

As Rouben has suggested, the use of Maple's VectorCalculus() package may be "overkill" for your requirements. A lot depends on just how much vector manipulation you plan on doing - which obvioulsy I don't know!. The attached shows simple addition of "free" vectors using the VectorCalculus() package

  restart;
  with(VectorCalculus):

  f:=(l, theta)->[ l*sin( convert(theta*degrees, radians)),
                   l*cos( convert(theta*degrees, radians))
                 ]:
  V1:= Vector( f(1, 30)):
  V2:= Vector( f(3/2, 60)):
  V3:=V1+V2;
  About(V3);
  PlotVector( [V1, V2, V3],
              color=[red, blue, green],
              width=[0.03, relative=false],
              head_length=[0.2, relative=false],
              axes=none
            );

Vector(2, {(1) = 1/2+(3/4)*sqrt(3), (2) = (1/2)*sqrt(3)+3/4})

 

Matrix(%id = 18446744074401223494)

 

 

 

 

 

 


 

Download simpleVects2.mw

@samiyare 

First of all, I probably should have spotted hat you were multiplying all values in the variable 't_number' by 60 - but I missed this, sorry. The effect of this "scaling" of the time axis, appears to just scale everything else - so one has to adjust ranges over whihc plots are made, search range for the optimizer, etc, etc. Note that both adjsuting ranges and controlling the optimizer range are used essentially to speed-up execution time: if you were pepared to wait quite a long time, then such range restrictions *might* be unnecessary - the relevant minimum might still be obtained, it would just take a loooonnngg time.

For what it is worth, the attached shows a minor variation of my original code implementing the 60*t_number with corresponding changes in plot ranges and optimization search space. This still locates the new "answers" in reasonable execution time


 

  restart;

#
# Parameters and procedures
#
  t_number:=<0, 10, 20, 30, 40>:
  m_number:=<11.50000000, 4.641182547, 1.273311993, 0.361845238, 0.288711649>:
  L:=2:
  getSol:= proc(d)
                local Mx0:= 0.05, cx0:= Mx0/(1-Mx0),
                      Mdb_i:= m_number[1], ct0:= Mdb_i,
                      PDE:= diff(C(x,t),t)=d*diff(C(x,t),x,x),
                      IBC:= {C(x,0)=ct0, C(0,t)=cx0, D[1](C)(L,t)=0};
                return pdsolve(PDE, IBC, numeric);
          end proc :
  Q:= proc(d)
           local pds, i, S:=0:
           if   not type(d, numeric)
           then return 'procname(_passed)'
           end if;
           pds:=getSol(d);
           for i from 1 to 5 do
               pds:-value( t=t_number[i]*60,
                           output=listprocedure
                         );                  # get the solution at the desired time
               eval(C(x,t), %);              # extract the C(x,t) at that time
               int(%, 0..2, numeric) / L;    # compute the average of C(x,t) at that time
               S += (% - m_number[i])^2;     # accumulate the residuals
           end do;
           return sqrt(S);
      end proc:

#
# Rouben's calculation
#
# A few representative values - shows that the
# "optimum" value for the parameter 'd' is
# somewhere around 0.14
#
# Change the following by ~10x to reflect the
# scaling factor
#
#  Q(0.1), Q(0.2), Q(0.5), Q(1.0);
#
  Q(0.001);Q(0.002);Q(0.003);Q(0.004);
  plot(Q(d), d=0.001 .. 0.004, numpoints=5, view=0..2);

HFloat(4.979783246316887)

 

HFloat(1.1695635853944253)

 

HFloat(1.5511287221444159)

 

HFloat(2.6009866748705806)

 

 

#
# Use the Optimization() package to get the value
# of diffusion coefficient which produces the
# minimum overall residual. Compute the solution
# at this "optimal" diffusion coefficient
#
  oSol:= Optimization:-Minimize(Q, 0..0.005);
  pdsOpt:= getSol( oSol[2][1] ):
#
# Using the optimal diffusion coefficient compute
# the average residual at each of the time-values
# in the vector t_number. These values *ought* to
# match somewhat with the corresponding value in
# the vector m_number
#
  resid:= Vector( 5,
                  i-> int
                      ( eval
                        ( C(x,t),
                          pdsOpt:-value
                                  ( t=t_number[i],
                                    output=listprocedure
                                  )
                        ),
                        0..2,
                        numeric
                      )/L
                 );
#
# Get the mean square error in the average C(x,t)-value
# at each of the time-values in the vector t_number
#
  `^`~((resid-m_number),2);

[.8889441947641348, Vector(1, {(1) = 0.2306692069320807e-2})]

 

Vector(5, {(1) = 11.499999999999842, (2) = 10.51283349395565, (3) = 10.108156818332095, (4) = 9.797236955370293, (5) = 9.534923705008406})

 

Vector[column](%id = 18446744074566322166)

(1)

#
# Plot C(x,t) for the optimal value of the diffusion
# coefficient for each of the time-values in t_number
#
  cols:= [red, green, blue, brown, black]:
  plots:-display( [ seq
                    ( pdsOpt:-plot
                      ( t=t_number[j],
                        x=0..L,
                        color=cols[j]
                      ),
                      j=1..numelems(t_number)
                    )
                  ]
                );
#
# Plot the "general 3-D" solution
#
  pdsOpt:-plot3d( x=0..L, t=0..t_number[-1]);

 

 

 


 

Download diffCoeff2.mw

you will have to upload the offending worksheet using the big green up-arrow in the Mapleprimes toolbar.

The TestDimensions() command i used by a lot of commands in the Units() system to verify that the units associated with terms in an expression are "compatible" - in other words making sure that on doesn't add "seconds" to "kilograms". A quick check shows that  it is capable of checking differetn ways of expressing the "same" unit, as shown in the attached "toy" example

restart;
kernelopts(version);
with(Units):
TestDimensions(2.114163508*10^7*Unit(kg/s^2)+2.114163508*10^7*Unit(J/m^2));

`Maple 2019.1, X86 64 WINDOWS, May 21 2019, Build ID 1399874`

 

Automatically loading the Units[Simple] subpackage

 

true

(1)

 


 

Download testUnit.mw

@mehran rajabi 

to the problem I answered in your other post here

https://www.mapleprimes.com/questions/227667-Solve-PDE-With-RungeKutta-Method-And-Plot

Posting the same question multiple times just confuses everyone - don't do it

@maple2015 

if anything- is the problem????

Your PDE is fourth-order in the variable 'x', and therefore needs 4 'x-related' boundary conditions - you have five

Your PDE is second-order in the variable 't', and therefore needs 2 't-related' initial conditions - you have one

The right hand side of this PDE (the piecewise function) looks a bit weird - it has a value at x=L, and absoutely nowhere else. This is probably going to result in a serious discontinuity an any solution (and may make the solution impossible to generate). Is this piecewise definition correct???

@Carl Love 

I have been reading this thread with increasing bewilderment. This confusion exists on two levels

  1. The worksheet  RLC.mw you posted does not generate the displayed answers! When I re-execute it (see attached below) the values of 'Sol' and related plot are different: note the arguments of the sin() and cos() terms should be 3*t, not 7*t. How is this possible? It is relatively easy to demonstrate that this coefficient (aka the the damped frequency) ought to be 3, although this demonstration is further confused by item 2 below
  2. As written, the ODE 'RLC' is incredibly confusing. You are using the symbol 'R' to denote inductance and the symbol 'L' to denote resistance, which is perverse to say the least. However it is the only way that applying Kirchoff's current law to to a series circuit of this type would make any sense. Of course one can use any names one wants, but OP's original post has R = 20 , C = 0.01 , L = 10, whereas  has   R=10, L=20, C= 1/100 so names have been swapped in the equation RLC and values have been swapped back in the params list. Hence the actual effect is nil (other than to confuse old people like me)
  3. In the second attachment below, I have defined the series RLC circuit in the conventional way, and swapped the parameter values around in order to achieve some kind of consistency. This means that both of these worksheets will now generate the same answers (but neither will generate the damped frequency value referred to in (1) above!)

restart:

RLC:= R*diff(x(t),t$2) + L*diff(x(t),t) + x(t)/C = 0;

R*(diff(diff(x(t), t), t))+L*(diff(x(t), t))+x(t)/C = 0

(1)

IC:= x(0)=10, D(x)(0)=0;

x(0) = 10, (D(x))(0) = 0

(2)

params:= [R=10, L=20, C= 1/100]:

Sol:= dsolve(eval({RLC, IC}, params));

x(t) = (10/3)*exp(-t)*sin(3*t)+10*exp(-t)*cos(3*t)

(3)

plot(eval(x(t), Sol), t= 0..5);

 

#Original square wave function:
step:= t-> piecewise(t < 0, 0, 1):
f:= t-> sum((-1)^n*step(t-n), n= 0..50):

#
#It's equivalent to this simpler function:
f:= t-> piecewise(floor(t)::even, 1, 0):

plot(f, 0..5, axes= frame, title= "Square Wave",thickness= 4, gridlines= false);

 

RLC2:= lhs(RLC) = 200*f(t/3);

RLC2 := R*(diff(x(t), t, t))+L*(diff(x(t), t))+x(t)/C = 200*piecewise((floor((1/3)*t))::even, 1, 0)

(4)

Sol2:= dsolve(eval({RLC2, IC}, params), numeric):

plots:-odeplot(Sol2, [[x(t), diff(x(t),t)]], t= 0..25, numpoints= 1000, thickness= 2, title= "Phase portrait", gridlines= false);

 

plots:-odeplot(Sol2, [t, x(t)], t= 0..25, numpoints= 1000, thickness= 2, title= "Solution curve", gridlines= false);

 

 

Download carlsCode.mw

restart:

#
# swapped round R and L
#
RLC:= L*diff(x(t),t$2) + R*diff(x(t),t) + x(t)/C = 0;

L*(diff(diff(x(t), t), t))+R*(diff(x(t), t))+x(t)/C = 0

(1)

IC:= x(0)=10, D(x)(0)=0;

x(0) = 10, (D(x))(0) = 0

(2)

#
# swapped rond values for R and L
#
params:= [R=20, L=10, C= 1/100]:

Sol:= dsolve(eval({RLC, IC}, params));

x(t) = (10/3)*exp(-t)*sin(3*t)+10*exp(-t)*cos(3*t)

(3)

plot(eval(x(t), Sol), t= 0..5);

 

#Original square wave function:
step:= t-> piecewise(t < 0, 0, 1):
f:= t-> sum((-1)^n*step(t-n), n= 0..50):

#
#It's equivalent to this simpler function:
f:= t-> piecewise(floor(t)::even, 1, 0):

plot(f, 0..5, axes= frame, title= "Square Wave",thickness= 4, gridlines= false);

 

RLC2:= lhs(RLC) = 200*f(t/3);

RLC2 := L*(diff(x(t), t, t))+R*(diff(x(t), t))+x(t)/C = 200*piecewise((floor((1/3)*t))::even, 1, 0)

(4)

Sol2:= dsolve(eval({RLC2, IC}, params), numeric):

plots:-odeplot(Sol2, [[x(t), diff(x(t),t)]], t= 0..25, numpoints= 1000, thickness= 2, title= "Phase portrait", gridlines= false);

 

plots:-odeplot(Sol2, [t, x(t)], t= 0..25, numpoints= 1000, thickness= 2, title= "Solution curve", gridlines= false);

 

 


 

Download convRLC.mw

@nm 

because we "froze" two different things. In your code

r    := (x^2-1)/x^2:
Z    := freeze(r):

whereas in mine

Z:=freeze( (x^2-1) ):

 

 

 

@nm 

one of the following more "suitable".

A "trick" I have observed with algsubs() is to (sometimes) do the substitution in a couple of stages


 

  restart;
  expr1:=(x^2-1)*y/x^2+(x^2-1)/x^2;

(x^2-1)*y/x^2+(x^2-1)/x^2

(1)

#
# The 'easy' way
#
  thaw(algsubs((y+1)=freeze((y+1)), expr1));

(y+1)*(x^2-1)/x^2

(2)

#
# Slightly harder - using algsubs and freeze/thaw
#
  Z:=freeze( (x^2-1) ):
  simplify
  ( collect
    ( algsubs
      ( x^2-1=Z,
        expr1
      ),
      Z
    )
  ):
  thaw(%);

(y+1)*(x^2-1)/x^2

(3)

#
# Even harder - using algsubs witout freeze/thaw
#
  subs
  ( Z1=(x^2-1),
    subs
    ( Z2=Z1/x^2,
      collect
      ( algsubs
        ( Z1/x^2=Z2,
          algsubs
          ( (x^2-1)=Z1,
            expr1
          )
        ),
        Z2
      )
    )
  );

(y+1)*(x^2-1)/x^2

(4)

 

 


 

Download frth3.mw

@minhthien2016 

to my previous worksheet to return the radius of the circumcircle and the incircle

triFun2.mw

odetest() only makes very careful simplifications when checking the validity of solutions.

It is sometimes(?) informative to perform the equivalent of odetest() manually - just substitute the solution into the ode, and make your own "simplifications", although I would advise doing this very carefully

For your specific example, the first command in the following evaluates to 'true', while the second evaluates to 'false' - the 'symbolic' option in the simplify() command makes all the difference

evalb(simplify(subs( my_sol, ode), symbolic));
evalb(simplify(subs( my_sol, ode)));

 

2 3 4 5 6 7 8 Last Page 4 of 147