tomleslie

9360 Reputation

17 Badges

11 years, 246 days

MaplePrimes Activity


These are answers submitted by tomleslie

You state

Maple code available to solve the heat  two dimensional heat equations using Rk4 

Invoking RK4 mean that you want to solve this numerically rather than analytically - correct? However as the help page for pdsolve/numeric states (my emphasis)

The pdsolve(PDEsys,conditions,numeric,vars,options) command returns a module that can be used to compute numerical solutions for time-based PDE systems over a fixed finite 1-space interval.

So a numerical solation over 2D-space is not (generally?) possible

It is however possible to generate an analytic solution over 2D space (see the attached). How explicit this can be made depends very much on the initial/boundary conditins of your problem, since it depends on 5 initial/boundary condiotns and 2 variable separation constants

restart;
PDE1 := diff(u(x, y, t), t)-diff(u(x, y, t), x$2)-diff(u(x, y, t), y$2) = 0;
Sol := pdsolve(PDE1, explicit);

diff(u(x, y, t), t)-(diff(diff(u(x, y, t), x), x))-(diff(diff(u(x, y, t), y), y)) = 0

 

u(x, y, t) = exp(_c[1]^(1/2)*x)*exp(_c[2]^(1/2)*y)*_C5*exp(t*_c[1])*exp(t*_c[2])*_C1*_C3+exp(_c[2]^(1/2)*y)*_C5*exp(t*_c[1])*exp(t*_c[2])*_C2*_C3/exp(_c[1]^(1/2)*x)+exp(_c[1]^(1/2)*x)*_C5*exp(t*_c[1])*exp(t*_c[2])*_C1*_C4/exp(_c[2]^(1/2)*y)+_C5*exp(t*_c[1])*exp(t*_c[2])*_C2*_C4/(exp(_c[1]^(1/2)*x)*exp(_c[2]^(1/2)*y))

(1)

 

Download heat.mw

 

ie

restart;
solve( { 0 = a_2*cos(7.5*b_2) + a_2, 250 = Pi*int((a_2*cos(b_2*x) + a_2)^2, x = 0 .. 7.5),
         750 = Pi*int((a_1*x^2 + b_1*x + sqrt(d_1*x) + 7.5)^2 - (a_2*cos(b_2*x) + a_2)^2, x = 0 .. x_2),
         2.5 = a_1*20^2 + 20*b_1 + sqrt(20*d_1) + 7.5,
         2*a_1*x_2 + b_1 + sqrt(d_1)/(2*sqrt(x_2)) = 0,
         (2*a_1)*20 + b_1 + sqrt(d_1)/(2*sqrt(20)) = 0
      },
      {a_1, a_2, b_1, b_2, d_1, x_2}
    )

It should be fairly obvious that

  1. the first three equations in the above do not contain x_2
  2. the last wo equations are identical - when x_2=20

The error message occurs because the integrand evaluates to a complex numbers for certain values of the integration variable 'x'. Some integration methods will handle complex integrands, but some won't

I'm not sure how meaningful it is, but if no integration method is specified, then Maple returns

A3 := Float(infinity) + 0.2486067977*I

see the attached

restart;
W := -16*(x - 1/2)*x*(t + 2*w)*(ln((m^2 + w*x*(-1 + x))/m^2)*ln(-1/(w*(t + w)*x^2 - w*(t + w)*x + m^2*t)) + ln(-x*w^2*(-1 + x))*ln((x^2 - x)*w + m^2) + ln(1/m^2)*ln((t + w)*(-1 + x)) + ln(-1/(x*w))*ln(m^2) + dilog(x*w^2*(-1 + x)/(w*(t + w)*x^2 - w*(t + w)*x + m^2*t)) - dilog(x*w*(-1 + x)*(t + w)/(w*(t + w)*x^2 - w*(t + w)*x + m^2*t)))/t^3 + 8*(x - 1/2)*((t + w)*x - t/2)*(-t*(w*x^2 + m^2 - w*x)*(t + w)*ln((x^2 - x)*w + m^2) + w*(x*w*(-1 + x)*(t + w)*ln((t + w)/w) + m^2*t*ln(m^2)))*(2*w*x + t)/(w*(t + w)*(w*(t + w)*x^2 - w*(t + w)*x + m^2*t)*t^3*x) + (16*x^2 - 8*x)/(t*w*(-1 + x)) + (2*t*x + 2*w*x - t - 2*w)*(2*w*x + t - 2*w)*(1 - 2*x)*ln((t*x + (-1 + x)*w)/((-1 + x)*w))/(((-1 + x)^2*w^2 + w*t*x*(-1 + x) + m^2*t)*t^2) + (-1 + 2*x)*(2*m^2 + w*x - w)*(2*w*x^2 + 2*m^2 - 3*w*x + w)*ln(m^2/(m^2 + w*x*(-1 + x)))/(w^2*(-1 + x)^2*((-1 + x)^2*w^2 + w*t*x*(-1 + x) + m^2*t)) + (-1 + 2*x)*(2*m^2 + w*x - w)*(2*w*x^2 + 2*m^2 - 3*w*x + w)*ln(m^2/(m^2 + w*x*(-1 + x)))/(w^2*(-1 + x)^2*(w*(t + w)*x^2 - w*(t + w)*x + m^2*t)) + (2*w*x + t)*((t + w)*x - t/2)*(x - 1/2)*ln(w^4/(t + w)^4)/(t^2*(w*(t + w)*x^2 - w*(t + w)*x + m^2*t)) - 8*(x - 1/2)*((-2 + 2*x)*w + t)*((-1 + x)*w + (x - 1/2)*t)*((t*x + (-1 + x)*w)*w^2*(-1 + x)^2*ln((t*x + (-1 + x)*w)/((-1 + x)*w)) + (-(t*x + (-1 + x)*w)*((x^2 - x)*w + m^2)*ln((x^2 - x)*w + m^2) + ln(m^2)*m^2*w*(-1 + x))*t)/((t*x + (-1 + x)*w)*w*(-1 + x)*t^3*((-1 + x)^2*w^2 + w*t*x*(-1 + x) + m^2*t)) + 16*(x - 1/2)*(ln(1/m^2)*ln((t*x + (-1 + x)*w)*(-1 + x)*w/((-1 + x)^2*w^2 + w*t*x*(-1 + x) + m^2*t)) + ln((x^2 - x)*w + m^2)*ln(1/((-1 + x)^2*w^2 + w*t*x*(-1 + x) + m^2*t)) + ln((-1 + x)^2*w^2)*ln((x^2 - x)*w + m^2) - dilog((t*x + (-1 + x)*w)*(-1 + x)*w/((-1 + x)^2*w^2 + w*t*x*(-1 + x) + m^2*t)) + dilog((-1 + x)^2*w^2/((-1 + x)^2*w^2 + w*t*x*(-1 + x) + m^2*t)))*(t*x + 2*(-1 + x)*w)/t^3:

f3 := (mm, ss, tt, ww) -> Int(eval(W, [s = ss, t = tt, w = ww, m = mm]), x = 0 .. 1):
A3 := evalf(f3(1, 3, -2, -1));

Float(infinity)+.2486067977*I

(1)

 


Download compInt.mw

 

for the basic equation in kinetics relating initial velocity, acceleration, time, and final velociy - in other words v=u+a*t, the simplest way is probably that shown in the attached

  restart;
  with(Units):
#
# The most basic equation in kinetics is probably
#
#     v=u*a*t
#
# So the initial velocity 'u' with a unit is
#
  initVel:=u*Unit(m/s);
#
# The acceleration with a unit is
#
  accel:=a*Unit(m/s/s);
#
# the time for whihc the acceleration is applied
# with  unit is
#
  tim:= t*Unit(s);
#
# So what is the final velocity (with units)
#
  finalVel:=initVel+accel*tim;
#
# And with numeric values......
#
  eval( finalVel, [u=2, a=3, t=4]);

Automatically loading the Units[Simple] subpackage
 

 

u*Units:-Unit(m/s)

 

a*Units:-Unit(m/s^2)

 

t*Units:-Unit(s)

 

(a*t+u)*Units:-Unit(m/s)

 

14*Units:-Unit(m/s)

(1)

 


 

Download simpU.mw

 

you can use the showstat() command.

However Maple commands come in two "flavours", those implemented in the kernel (also known as "builtin" commands) and those in the "library". You cannot view the code for any "builtin" command, but showstat() will display the code for library commands.

Unfortunately both taylor and laurent are "builtin" commands, so attempting showstat(taylor) or showstat(laurent) will hust return the error

Error, (in showstat) cannot debug built-in functions

You can contrast this with (for example) showstat(mtaylor) which will display the code implemented by the command.

Note that the code returned by showstat() may (probably will!) contain references of "builtin" commands

depending on where the data is coming from, any restrictions on the data etc, etc

One possibility for a "toy" example is shown in the attached.

Not sure why the resulting figure does not rendeer properly on this site - it does in the worksheet -honest!

  restart;
  with(plottools):
  with(plots):
  r:=rand(1..20):
  xvals:=[$i=1..20]:
  yvals:=[seq( r(), j=1..20)]:
  rects:=[seq( rectangle( [xvals[j]-1/2, 0], [xvals[j]+1/2, yvals[j]], color=blue), j=1..20)]:
  display(rects);

 

 

Download histo.mw

 

Your original worksheet appears to do more or less the same calculation multiple times (maybe with somewhat different parameters?)

The attached does only the first "loop" calculation and produces the attached Excel file.

You will have to change the filepath in the final ExcelTools:-Export() command to something appropriate for your installation

aCalc.mw

resMat.xlsx

  1. It is your responsibility to ensure that matrix dimensions are "compatible" for which ever matrix operation you are attempting to perform.
  2. There can be an "issue" over the interpretation of the '.' operation if you mix Vector() and Matrix() definitions. If you are performing "conventional" matrix operations, then it is probably safer to define everything as a Matrix(), even for 1xn or nx1 structures. If you use Vector() definitions then it is possible under certain circumstances that the '.' operator will be interpreted as a dot (ie scalar) product, and hence produce a scalar quantity, rather than a 1x1 matrix

Contemplate the difference between V1.V2 and V3.V4 in the attached!

  restart:
  V1:= Matrix(1, 3, [1, 2, 3]):
  V2:= Matrix(3, 1, [4, 5, 6]):
  V1.V2; # outputs a 1x1 Matrix
  V2.V1; # outputs a 3x3 Matrix
  V3:= Vector[row]([1, 2, 3]):
  V4:= Vector[column]([4, 5, 6]):
  V3.V4;  # performs the scalar (ie "dot") product so outputs a scalar
  V4.V3;  # outputs a 3x3 Matrix

Vector(1, {(1) = 32})

 

Matrix(3, 3, {(1, 1) = 4, (1, 2) = 8, (1, 3) = 12, (2, 1) = 5, (2, 2) = 10, (2, 3) = 15, (3, 1) = 6, (3, 2) = 12, (3, 3) = 18})

 

32

 

Matrix(%id = 36893488147955157220)

(1)

 

Download matVec.mw

 

 

see the attached

  restart;
  f := theta -> sin(x(theta));
  D(f)(t);
  convert(%, diff);

proc (theta) options operator, arrow; sin(x(theta)) end proc

 

(D(x))(t)*cos(x(t))

 

(diff(x(t), t))*cos(x(t))

(1)

  restart;
  f:=sin(x(t));
  diff(f, t);

sin(x(t))

 

(diff(x(t), t))*cos(x(t))

(2)

 

Download diffs.mw

try using 100 raised to the power o (as in the letter 'o') ie 100^o

What do you want to happen when more than one element satisfied the criterion??

While you think about that, you might consider the "toy" example shown in the attached whch shows a couple of different ways

  restart;
  L:=[seq( j/2, j=10..0,-1)];
  seq( `if`( `and`(L[j]>0, L[j]<1), [j,L[j]], NULL), j=1..numelems(L));
  ListTools:-SelectFirst(numelems(L), j-> `and`(j>0, j<1), L, output=[indices, values]);

[5, 9/2, 4, 7/2, 3, 5/2, 2, 3/2, 1, 1/2, 0]

 

[10, 1/2]

 

[10], [1/2]

(1)

 

Download listSel.mw

This is known as the "Context Bar"

So make sure it is toggled in the View menu

is shown in the attached. (The plot renders rather differently in the worksheet thsn it does on this site!)

  restart;
  with(Statistics):
  x:= RandomVariable(Binomial(45, 0.9)):
  data:=[seq(ProbabilityFunction(x, i), i=0..45)]:
  dataplot(data);

 

 

Download plFata.mw

When you compute the "same" quantity in two completely different ways, you really shouldn't be surprised if the answers "look" different. The important question is - are they different?

The attached shows that both methods  in fact give the same answer, after appropriate simplification. Check the final entry from each execution group.

  restart:
  with(Syrup):
  with(inttrans):
#
# Set up the problem
#
  ckt:= [ vin(5), Rm(1.00000000*10^7), Cm(1.80000000*10^(-11), ic = 0),
          R1(9.00000000*10^6) &// C1(1.50000000*10^(-11), ic = 0),
          R2(1.00000000*10^6) &// C2(1.50000000*10^(-11), ic = 0)
        ]:
  Draw(ckt);

 

#
# Get the voltage across the capacitor by solving
# the odes (basically a transient analysis)
#
  vcm:= eval
        ( v[Cm](t),
          dsolve
          ( Solve
            ( ckt,
              'tran'
            )
          )
        );
  simplify(expand(vcm));

  

(1/3)*exp((100000/459)*(-119+8041^(1/2))*t)*((61/64328)*8041^(1/2)-1/8)*8041^(1/2)-(1/3)*exp(-(100000/459)*(119+8041^(1/2))*t)*(-(61/64328)*8041^(1/2)-1/8)*8041^(1/2)+(91/3)*exp((100000/459)*(-119+8041^(1/2))*t)*((61/64328)*8041^(1/2)-1/8)+(91/3)*exp(-(100000/459)*(119+8041^(1/2))*t)*(-(61/64328)*8041^(1/2)-1/8)+5/2

 

(1/32164)*(-415*8041^(1/2)-40205)*exp((100000/459)*(-119+8041^(1/2))*t)+5/2+(1/32164)*(415*8041^(1/2)-40205)*exp(-(100000/459)*(119+8041^(1/2))*t)

(1)

#
# Get the transfer function to the node '2' (basically
# an ac analysis) then compute the step response, and
# convert the result from hyperbolic trig functions to
# exponentials
#
  v2:= convert
       ( invlaplace
         ( eval
           ( v[2],
             [ Solve
               ( ckt,
                 returnall
               )
             ][1]
           )/s,
           s,
           t
        ),
        exp
      );
  simplify(expand(v2));

5/2-(5/32164)*(83*exp((200000/459)*t*8041^(1/2))*8041^(1/2)+8041*exp((200000/459)*t*8041^(1/2))-83*8041^(1/2)+8041)*exp(-(700000/27)*t-(100000/459)*t*8041^(1/2))

 

(1/32164)*(-415*8041^(1/2)-40205)*exp((100000/459)*(-119+8041^(1/2))*t)+5/2+(1/32164)*(415*8041^(1/2)-40205)*exp(-(100000/459)*(119+8041^(1/2))*t)

(2)

 

Download ckt.mw

of which one is shown in the attached

  restart;
  with(RootFinding):
  sols:= []:
  tmp:=0:
  while true do
        tmp:= NextZero
              (  x -> sin(x)-cos(4*x - 1/6*Pi),
                 tmp
              );
        if   tmp <= evalf(2*Pi) then
             sols:= [sols[], tmp];
        else break;
        fi:
  end do:
  sols;

[.4188790204, 1.675516081, 1.745329251, 2.932153143, 3.839724354, 4.188790204, 5.445427266, 5.934119456]

(1)

 

Download roots.mw

 

1 2 3 4 5 6 7 Last Page 1 of 166