tomleslie

Hmm...

Answered: tomleslie 13876

November 01 2021

1 1

I can solve the reverse problem, with a couple of extra assumptions - this may (or may not!) be of some help

>	restart; f := beta/(2sigmaGAMMA(1/beta))*exp(-(abs(x-mu)/sigma)^beta); assume(mu::real, beta > 0, sigma > 0); F := int(f, x=-infinity..s);

(1/2)*beta*exp(-(abs(-x+mu)/sigma)^beta)/(sigma*GAMMA(1/beta))

int((1/2)*beta*exp(-(abs(-x+mu)/sigma)^beta)/(sigma*GAMMA(1/beta)), x = -infinity .. s)

(1)

>	assume(mu>0, x>0): # unjustified assumptions F := 1/2+signum(x-mu)/(2GAMMA(1/beta))(GAMMA(1/beta)-GAMMA(1/beta, abs((x-mu)/sigma)^beta)); diff(F, x); ans:=simplify(%, symbolic); is(ans=f);

1/2-(1/2)*signum(-x~~+mu~~~)*(GAMMA(1/beta)-GAMMA(1/beta, (abs(-x~~+mu~~~)/sigma)^beta))/GAMMA(1/beta)

-(1/2)*signum(1, -x~~+mu~~~)*(GAMMA(1/beta)-GAMMA(1/beta, (abs(-x~~+mu~~~)/sigma)^beta))/GAMMA(1/beta)+(1/2)*signum(-x~~+mu~~~)^2*(abs(-x~~+mu~~~)/sigma)^beta*beta*((abs(-x~~+mu~~~)/sigma)^beta)^(1/beta-1)*exp(-(abs(-x~~+mu~~~)/sigma)^beta)/(GAMMA(1/beta)*abs(-x~~+mu~~~))

(1/2)*beta*exp(-abs(-x~~+mu~~~)^beta*sigma^(-beta))/(sigma*GAMMA(1/beta))

(2)

>

Download manip.mw

A few things...

Answered: tomleslie 13876

November 01 2021

2 9

in the definition of your pde, you use psi(x), but 'psi' is an expression, not a function - so I have corrected this.
in the definition of psi, two conditions apply when x=l/2, since both evaluate to the same thing this probably isn't too important, but you really shouldn't do this.
when solving the wave equation it is always a good idea to provide the pdsolve() command with the 'hint=TWS' option
I have inserted the assumption that l>0. Without this assumption, Maple provides three solutions, one for l<0, one for l=0, and one for l>0. The first two of these are not very interesting/useful

The attached returns an explicit solution in the form which you want (and it is simpler!)

(I'm not sure why this website is displaying Maple output in such an untidy way! I think something on this site is "broken")

>

  restart;
  assume( l>0 ):
  interface(showassumed=0):
# set the PDE, ic(initial conds) and bc (boundary conds)

  pde := diff(u(x, t), t, t) = a^2*diff(u(x, t), x, x);
  ic := u(x, 0) = 0, D[2](u)(x, 0) = psi;
  bc := u(0, t) = 0, u(l, t) = 0;
#
# OP's original definition of psi had a typo
#
  psi := piecewise( 0 <= x and x <= l/2, x,
                    l/2 < x and x <= l, l - x
                  );

  res := pdsolve([pde,ic,bc], HINT = TWS);

piecewise(0 <= x and x <= (1/2)*l, x, (1/2)*l < x and x <= l, l-x)

u(x, t) = Sum(4*l^2*sin((1/2)*n*Pi)*sin(n*Pi*x/l)*sin(a*n*Pi*t/l)/(Pi^3*n^3*a), n = 1 .. infinity)

(1)

>

Download waveEq.mw

Three things...

Answered: tomleslie 13876

October 29 2021

3 6

The command is dsolve() not desolve()
I don't think dsolve() handles expressions with units so I converted you initial condition to be unit-ffree
No value is specified for P_b so I (arbitrarily) chose 1.0

and everything in the attached now seems to "work"

Automatically loading the Units[Simple] subpackage

Variables

Physical parameters

(1.1.1)

(1.1.2)

(1.1.3)

(1.1.4)

(1.1.5)

(1.1.6)

Variables

(1.2.1)

(1.2.2)

(1.2.3)

(1.2.4)

(1.2.5)

(1.2.6)

(1.2.7)

Calculate water vapor

The formula used below comes from Huang 2018 "A Simple Accurate Formula for Calculating Saturation Vapor Pressure of Water and Ice"

(2.1)

"`P__s`(T):=(exp(34.494-4924.99/((T+237.1))))/((T+105)^(1.57))"

(2.2)

(2.3)

(2.4)

Temperature conversion PO2

(3.1)

(3.2)

(3.3)

(3.4)

(3.5)

(3.6)

(3.7)

(3.8)

"S(PO2,P50):=(((PO2/(P50))^(n))/(1+(PO2/(P50))^(n)))"

(3.9)

(3.10)

(3.11)

(3.12)

(3.13)

"`PO2__c`(S,P50):=(exp((ln((-S)/((S-1))))/(n))*P50)"

(3.14)

(3.15)

(3.16)

Calculation Mass Transfer

Solubility

k__c := 4.658*10^(-5)*1.01^(37-T__b)*Unit('mL'/'mmHg'*(1/'mL'))

(4.1.1)

k__p := 2.855*10^(-5)*1.01^(37-T__b)*Unit('mL'/'mmHg'*(1/'mL'))

(4.1.2)

(4.1.3)

"lambda(`PvO2__c`):=(1.34/(convert(k,unit_free))*convert(Hgb,unit_free) *(n/(convert(`P50__c`,unit_free)))*((convert(`PvO2__c`,unit_free))/(convert(`P50__c`,unit_free)))^((n-1))*(1/((1+((convert(`PvO2__c`, unit_free))/(convert(`P50__c`,unit_free)))^(n))^(2))))/(1000)"

(4.1.4)

(4.1.5)

Diffusivity

(4.2.1)

(4.2.2)

N := k__c/k__p

(4.2.3)

beta := (1/3)*(2/(1+(N*D__c/D__p-1)*(.283*(1/2)))+1/(1+(N*D__c/D__p-1)*(1-.283)))*(N*D__c/D__p-1)

(4.2.4)

`ϰ` := -(N*D__c/D__p-1-N*D__c*beta/D__p)/(N*D__c/D__p-1-beta)

(4.2.5)

R := Hct*(N*D__c/D__p-1)/(N*D__c/D__p+`ϰ`)

(4.2.6)

(4.2.7)

Deff := D__m/(1+lambda(PvO2__c))

(4.2.8)

Density

(4.3.1)

(4.3.2)

(4.3.3)

Viscosity

(4.4.1)

(4.4.2)

(4.4.3)

(4.4.4)

Dimensionless Numbers

N__re := Q__blood*d/((1-Void__f)*A__f*nu)

(4.5.1)

(4.5.2)

Differential Equation

Convert flow from mL/min to mL/s

(4.6.1)

Convert η from mPa*s to g/cm*s

(4.6.2)

a := 4*((1-Void__f)/convert(d, unit_free))^(1+`ϕ`)*(convert(A__f, unit_free)*convert(`η__c`, unit_free)/(convert(Q__bc, unit_free)*convert(rho, unit_free)))^`ϕ`*phi/(Void__f*(nu/D__m)^(2/3))

(4.6.3)

ODE := diff(P(x), x) = a*(P__b-P(x))/(1+lambda(P(x))^(2/3))

diff(P(x), x) = .1081111001*(P__b-P(x))/(1+1.846187058*(P(x)^1.7/(1+0.2439243816e-3*P(x)^2.7)^2)^(2/3))

(4.6.4)

(4.6.5)

$sol := dsolve(eval({IC, ODE}, P__b = 1.0), numeric); plots:-odeplot(sol, [x, P(x)], x = 0 .. 100)$

Download odeProb.mw

Hmmm...

Answered: tomleslie 13876

October 25 2021

4 1

Lists are not mutable containers so it is actually very inefficient to add elements one-by-one. Each time you add a new element you are actually creating a whole new list. As a purely practical matter, you probably won't notice this inefficiency until your list get really long. But it is a bad habit to get into!

On the other hand Arrays are mutable and can be grown (in-place) as desired by using "round-bracket" indexing. If necessary the end result can be converted to a list. See the attached

Download arrList.mw

Many ways to do this...

Answered: tomleslie 13876

October 25 2021

3 0

of whihc the attached show only one.

This worksheet shows that

if m=3, the roots of the expression are [3, 9]
if m=-2, the roots of the expression are [-2, 4]

>

  restart;
#
# Assume the OP only wants real solutions
#
  with(RealDomain):
#
# Insert OP's expression and solve for 'x'
# in terms of 'm'
#
  expr:= x^2+2*m*x-6*x=m^3;
  sol:= [ solve(expr, x) ];
#
# Find all values of 'm' in the above solutions
# where one root is the square of the other root
#
  mvals:= { solve( sol[1]=sol[2]^2, m),
             solve( sol[2]=sol[1]^2, m)
           };
#
# Find the roots of the original expression for
# the values obtained for 'm' above
#
  rts1:=solve(eval(expr, m=mvals[1]), x);
  rts2:=solve(eval(expr, m=mvals[2]), x);

(1)

>

Download sqRoots.mw

Well...

Answered: tomleslie 13876

October 22 2021

3 2

I would have expected your code to "work", since (in general) a list of linestyles is an allowed plot option - however in the case of DEPlot(), as you have observed, this cuases an error. This may be a bug(ette).

One can "work around" the issue by creating each plot individually and then displaying them together, as in the attached.

On your second issue, the 'axes' option accepts a single value - so you are stuck with either "boxed" or "normal". It is possible using a combination of textplot() and line plots to create a "plot" which "looks" like axes. This could be included in the display() command in the attached.

Note that although "gridlines" appear on the plot when rendered on this site, they do not occur in the worksheet (honest!)

>

  restart;
  with(VectorCalculus):
  with(LinearAlgebra):
  with(DEtools):
  F0 := 0:
  zeta := .25:
  w := 1:
  Omega := 1:
  m := 1:

  ode1 := diff(x(t), t) = y(t):
  ode2 := diff(y(t), t) = F0*cos(Omega*t)/m-2*zeta*w*y(t)-w^2*x(t):

  inits:=[ [x(0) = 0, y(0) = 50],
           [x(0) = 9, y(0) = 25],
           [x(0) = 85, y(0) = 20],
           [x(0) = .25, y(0) = .5],
           [x(0) = 7, y(0) = 5]
         ]:
  ls:= [dash, dashdot, dot, longdash, spacedash]:
  cols:= [red, green, black, navy, maroon]:
  p1:= [ seq
         ( DEplot
            ( [ode1, ode2],
              [x(t), y(t)],
              t = 0 .. 10,
              x = -100 .. 100,
              y = -100 .. 100,
              [ inits[k] ],
              linestyle=ls[k],
              linecolor =cols[k],
              thickness = 5,
              axes = normal,
              labels = [x(t), y(t)],
              arrows = small,
              color=cyan,
              stepsize = 0.1e-1,
              axes=boxed
           ),
           k=1..5
        )
      ]:
  plots:-display(p1);

>

Download dplt.mw

Hmmm - some observations...

Answered: tomleslie 13876

October 18 2021

1 1

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

Download heat.mw

If I write your equaions "neatly"...

Answered: tomleslie 13876

October 18 2021

2 1

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

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

Well...

Answered: tomleslie 13876

October 17 2021

2 3

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

(1)

>

Download compInt.mw

Simpler?...

Answered: tomleslie 13876

October 15 2021

0 0

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

(1)

>

Download simpU.mw

Generally speaking...

Answered: tomleslie 13876

October 12 2021

1 1

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

Many ways to do this...

Answered: tomleslie 13876

October 12 2021

1 0

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!

Download histo.mw

One way is shown in the attached...

Answered: tomleslie 13876

October 12 2021

1 2

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

Two suggestions...

Answered: tomleslie 13876

October 11 2021

1 0

It is your responsibility to ensure that matrix dimensions are "compatible" for which ever matrix operation you are attempting to perform.
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!