C_R

3932 Reputation

21 Badges

6 years, 319 days

MaplePrimes Activity


These are questions asked by C_R

In the below, the example throws an error

restart;
foo := proc(params) 
local names; 
names := lhs~(params)[]; 
#local op(names);
print([names]):
assign(params); 
print(b, c, d); end proc;

foo([b = 2, c = 3, d = 4]);
foo([b = 2, c = 3, d = 4]); #  error because the assignement was done to global names

How to fix that and/or avoid assignements to global names

For a numericial integration I want to check if all physical parameters are given in SI base units (before stripping of the units for the integration). Right now I am doing this

restart

Parameters

H = 0.1e-1*Unit('m'), R = 0.6e-2*Unit('m'), r_sub = 0.3e-2*Unit('m'), k = .540*Unit('W'/('m'*'K')), rho = 1060.0*Unit('kg'/'m'^3), Cp = 3745*Unit('J'/('kg'*'K')), Q_flux = 5000*Unit('W'/'m'^2), T_amb = 0*Unit('K'), Nz = 20, Nr = 20, t_final = 30*Unit('s'), t_step = Unit('s'), test_1 = Unit('kW'), test_2 = 2*Unit('N')

H = 0.1e-1*Units:-Unit(m), R = 0.6e-2*Units:-Unit(m), r_sub = 0.3e-2*Units:-Unit(m), k = .540*Units:-Unit(W/(m*K)), rho = 1060.0*Units:-Unit(kg/m^3), Cp = 3745*Units:-Unit(J/(kg*K)), Q_flux = 5000*Units:-Unit(W/m^2), T_amb = 0, Nz = 20, Nr = 20, t_final = 30*Units:-Unit(s), t_step = Units:-Unit(s), test_1 = Units:-Unit(kW), test_2 = 2*Units:-Unit(N)

(1)

H: length of rod
R: outer radius
r_sub: outer radius of heat flux
k: conductance
ρ: density
Cp: heat capacity
Q_flux: heat flux
T_amb: initial Tempearture
Nz: number of axial grid points
Nr: number of radial grid points
t_final: integration time
t_step: sampling interval for temperature fields
test_1: unit given in kilos
test_2: unit given in derived base units

alpha = k/(rho*Cp); simplify(subs(H = 0.1e-1*Units:-Unit(m), R = 0.6e-2*Units:-Unit(m), r_sub = 0.3e-2*Units:-Unit(m), k = .540*Units:-Unit(W/(m*K)), rho = 1060.0*Units:-Unit(kg/m^3), Cp = 3745*Units:-Unit(J/(kg*K)), Q_flux = 5000*Units:-Unit(W/m^2), T_amb = 0, Nz = 20, Nr = 20, t_final = 30*Units:-Unit(s), t_step = Units:-Unit(s), test_1 = Units:-Unit(kW), test_2 = 2*Units:-Unit(N), %))

alpha = 0.1360304305e-6*Units:-Unit(m^2/s)

(2)

α: thermal diffusivity

 

One way is to check whether the value of a physical parameter has not changed after simplification

`minus`(map(convert, simplify(`~`[rhs]({Cp = 3745*Units:-Unit(J/(kg*K)), H = 0.1e-1*Units:-Unit(m), Nr = 20, Nz = 20, Q_flux = 5000*Units:-Unit(W/m^2), R = 0.6e-2*Units:-Unit(m), T_amb = 0, k = .540*Units:-Unit(W/(m*K)), r_sub = 0.3e-2*Units:-Unit(m), rho = 1060.0*Units:-Unit(kg/m^3), t_final = 30*Units:-Unit(s), t_step = Units:-Unit(s), test_1 = Units:-Unit(kW), test_2 = 2*Units:-Unit(N)})), unit_free), map(convert, `~`[rhs]({Cp = 3745*Units:-Unit(J/(kg*K)), H = 0.1e-1*Units:-Unit(m), Nr = 20, Nz = 20, Q_flux = 5000*Units:-Unit(W/m^2), R = 0.6e-2*Units:-Unit(m), T_amb = 0, k = .540*Units:-Unit(W/(m*K)), r_sub = 0.3e-2*Units:-Unit(m), rho = 1060.0*Units:-Unit(kg/m^3), t_final = 30*Units:-Unit(s), t_step = Units:-Unit(s), test_1 = Units:-Unit(kW), test_2 = 2*Units:-Unit(N)}), unit_free))

{1000}

(3)

if is(`minus`(map(convert, simplify(`~`[rhs]({Cp = 3745*Units:-Unit(J/(kg*K)), H = 0.1e-1*Units:-Unit(m), Nr = 20, Nz = 20, Q_flux = 5000*Units:-Unit(W/m^2), R = 0.6e-2*Units:-Unit(m), T_amb = 0, k = .540*Units:-Unit(W/(m*K)), r_sub = 0.3e-2*Units:-Unit(m), rho = 1060.0*Units:-Unit(kg/m^3), t_final = 30*Units:-Unit(s), t_step = Units:-Unit(s), test_1 = Units:-Unit(kW), test_2 = 2*Units:-Unit(N)})), unit_free), map(convert, `~`[rhs]({Cp = 3745*Units:-Unit(J/(kg*K)), H = 0.1e-1*Units:-Unit(m), Nr = 20, Nz = 20, Q_flux = 5000*Units:-Unit(W/m^2), R = 0.6e-2*Units:-Unit(m), T_amb = 0, k = .540*Units:-Unit(W/(m*K)), r_sub = 0.3e-2*Units:-Unit(m), rho = 1060.0*Units:-Unit(kg/m^3), t_final = 30*Units:-Unit(s), t_step = Units:-Unit(s), test_1 = Units:-Unit(kW), test_2 = 2*Units:-Unit(N)}), unit_free)) <> {}) then print("Not all units are in SI base units") end if

"Not all units are in SI base units"

(4)

Download check_for_SI_base_units.mw

I am interested in simpler ways to perform the test. The statement is quite long and not easy to understand.

Given the size of the unit package, there might be a command that I have overlooked to check whether units are in SI base units.

Have I missed test cases to make the check fool proof?

 

Update: A variant with split was added and attempts with functional programming

check_for_SI_base_units_02.mw

This piecewise function with units

q := r -> piecewise(0 <= r and r < 0.1*Unit('mm'), 5*Unit(('kW')/'m'^2));

returns this output

Is it possible to get the output printed without ticks and the units in roman?

In the context of numerical integration of partial differential equations I am looking for tools to create difference equations for the numerical integration of PDE's. Ideally the desired integration scheme could be given as an input.

The only command that I am aware of is the DiffEquation command from the DynamicSystems package (see ?DynamicSystems,DiffEquation) which is limited to linear time invariant systems and only provides a simple forward integration scheme.

Any insights/examples into the topic of numerical integration of pde's with Maple beyond pdsolve/numeric capabilites is very much appreciated. (I know that Maple comes with advanced numerical solvers for ODE's. MapleSim uses such solvers to integrate systems of ODE's. This seems to be, in terms of computation, close to what I am looking for.)

I am pretty new to pdsolve.

?pdsolve/numeric states

PDEsys - single or set or list of time-dependent partial differential equations in two independent variables

This could be interpreted in a way that problems with 3 independent variables can be handled (time + two others).

The rest of the helppage only describes two independent variables and does not mention time.

In case the above interpretation is wrong, I have two questions:

  • If only two independent variables can be treated, why is time-depencency mentioned at all?
  • Is it possible to integrate the PDE of an axially heated rod of finite length and radius with the help of pdsolve where on one end of the rod a heat flux q=q(r) is applied?

Any links to similar problems solved with Maple (not necessary with pdsolve) are very much appreciated.

Edit: Updated BCs

PDE := diff(T(r, z, t), t) = alpha*(diff(T(r, z, t), r, r)+(diff(T(r, z, t), r))/r+diff(T(r, z, t), z, z))

diff(T(r, z, t), t) = alpha*(diff(diff(T(r, z, t), r), r)+(diff(T(r, z, t), r))/r+diff(diff(T(r, z, t), z), z))

(1)

Updated BC

BC1 := -k*(D[2](T))(r, 0, t) = q; BC2 := (D[2](T))(r, L, t) = 0; BC3 := (D[1](T))(R, z, t) = 0; BC4 := (D[1](T))(0, z, t) = 0; IC := T(r, z, 0) = 0

T(r, z, 0) = 0

(2)

NULL

Corrected original BC (this is a different case of constant external temperature)

BC1 := -k*(D[2](T))(r, 0, t) = q; BC2 := T(r, L, t) = 0; BC3 := T(R, z, t) = 0; BC4 := (D[1](T))(0, z, t) = 0; IC := T(r, z, 0) = 0

T(r, z, 0) = 0

(3)

NULL

Download 2d_BVP.mw

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