jetboo

45 Reputation

3 Badges

2 years, 108 days

MaplePrimes Activity


These are replies submitted by jetboo

@mmcdara

d1 is a centered divided difference approximation of diff(u, x$3) at "point" [i+1/2, j];not an approximation of the laplacian of u at [i, j]

yes you are right, I am using a 2D finite volume code and I need to evaluate on the faces of given cell  grad(laplacian(u)).n*length_of_face where n is the outward normal. On the right face of a cell located at i+1/2,j I need to evaluate diff(u,x,x,x) and diff(u,y,y,x) my comment was about the first part only.

Could you please comment on the syntax you used for the definition of tau ? I do not understand the double right arrow and the *~ operation.

Thank you again

@tomleslie 

thank you very much !

Does anyone have an  idea on how to assign symbolically a series of values, say alpha_n, to be the (infinite number of) roots of a Bessel function ? 

Cheers

I am trying to solve the heat conduction problem in an axisymmetric cylindrical bar with Laplace transforms, for which an analytical solution exists. I end up finding the right solution in Laplace space but I can't compute its inverse as I need to use the inversion theorem and compute residuals at poles. It implies that I have to assign a series of values, say alpha_n to be the roots of the Bessel function (of the first kind or order 0) J_0.  Is there a way to do that ?  Here is my maple page that shows the problem 

restart; with(PDEtools)

NULLwith(inttrans)

assume(constant*K_1 > 0)

 

about(K_1)

Originally K_1, renamed K_1~:
  Involved in the following expressions with properties
    constant*K_1 assumed RealRange(Open(0),infinity)
  also used in the following assumed objects
  [constant*K_1] assumed RealRange(Open(0),infinity)
 

 

# Equations for c in physical space, we want to solve the heat equation in an axi-symmetric cylindrical barPDE1 := diff(c(r, t), t) = K_1*(diff(r*(diff(c(r, t), r)), r))/r

diff(c(r, t), t) = K_1*(diff(c(r, t), r)+r*(diff(diff(c(r, t), r), r)))/r

(1)

NULL

ic1 := c(r, 0) = 0

c(r, 0) = 0

(2)

  NULLNULL

PDEL1 := laplace((lhs-rhs)(PDE1), t, p)

p*laplace(c(r, t), t, p)-c(r, 0)-K_1*laplace(diff(diff(c(r, t), r), r), t, p)-K_1*laplace(diff(c(r, t), r), t, p)/r

(3)

``

# Apply initial condition to simplify the expression

simplified_PDEL1 := eval(PDEL1, ic1)

p*laplace(c(r, t), t, p)-K_1*laplace(diff(diff(c(r, t), r), r), t, p)-K_1*laplace(diff(c(r, t), r), t, p)/r

(4)

NULL
eq_C := subs(laplace(c(r, t), t, p) = C(r, p), laplace(diff(c(r, t), r, r), t, p) = diff(C(r, p), r, r), laplace(diff(c(r, t), r), t, p) = diff(C(r, p), r), simplified_PDEL1)

p*C(r, p)-K_1*(diff(diff(C(r, p), r), r))-K_1*(diff(C(r, p), r))/r

(5)

# Boundary condition at r = a

bc_right := C(a, p) = V/p

C(a, p) = V/p

(6)

# One can directly inject bc_right in the solver dsolve but we choose here not toSol := dsolve(eq_C, C(r, p))

C(r, p) = _F1(p)*BesselJ(0, (-p/K_1)^(1/2)*r)+_F2(p)*BesselY(0, (-p/K_1)^(1/2)*r)

(7)

NULL

``

NULL (Bessel of the second kind of order 0 tends to infinity as r->0)

bounded_solution_C := subs(coeff(rhs(Sol), BesselY(0, sqrt(-p/K_1)*r)) = 0, Sol)NULL

C(r, p) = _F1(p)*BesselJ(0, (-p/K_1)^(1/2)*r)

(8)

#Apply now boundary condition at r=a given above

simplify(subs(r = a, rhs(bounded_solution_C)) = V/p)

_F1(p)*BesselJ(0, (-p/K_1)^(1/2)*a) = V/p

(9)

isolate(_F1(p)*BesselJ(0, (-p/K_1)^(1/2)*a) = V/p, _F1(p))

_F1(p) = V/(p*BesselJ(0, (-p/K_1)^(1/2)*a))

(10)

bounded_solution_C := subs(_F1(p) = V/(p*BesselJ(0, (-p/K_1)^(1/2)*a)), bounded_solution_C)

C(r, p) = V*BesselJ(0, (-p/K_1)^(1/2)*r)/(p*BesselJ(0, (-p/K_1)^(1/2)*a))

(11)

is(BesselJ(0, sqrt(-p/K_1)*r) = BesselJ(0, sqrt(-p/K_1)*r))

true``

(12)

`assuming`([c(r, t) = invlaplace(rhs(simplify(bounded_solution_C)), p, t)], [t > 0, K_1 > 0])

c(r, t) = V*invlaplace(BesselJ(0, (-p/K_1)^(1/2)*r)/(p*BesselJ(0, (-p/K_1)^(1/2)*a)), p, t)

(13)

NULL and is thus incomplete

`` 

is(BesselI(0, I*x) = BesselJ(0, x))

true

(14)

is(BesselJ(0, I*x) = BesselI(0, x))

true

(15)

is(BesselJ(0, I*x) = BesselI(0, I*x))

false``

(16)

# Inversion theorem for the Laplace transforms

"c(r,t) := 1/((2 Pi i) )(∫)[ gamma - i infinity]^(gamma +i infinity)V( (e)^((lambda t)) I[0](sqrt(lambda/(K_1)) r ))/(lambda I[0](sqrt(lambda/(K_1)) a)) ⅆlambda"

proc (r, t) options operator, arrow, function_assign; -((1/2)*I)*(int(V*exp(lambda*t)*`#msub(mn("I",fontstyle = "italic"),mn("0"))`(sqrt(lambda/K_1)*r)/(lambda*`#msub(mn("I",fontstyle = "italic"),mn("0"))`(sqrt(lambda/K_1)*a)), lambda = gamma-infinity*I .. gamma+infinity*I))/Pi end proc

(17)

# Pole at lamba = 0

NULL

integrand := exp(lambda*t)*BesselI(0, sqrt(lambda/K_1)*r)/BesselI(0, sqrt(lambda/K_1)*a)

exp(lambda*t)*BesselI(0, (lambda/K_1)^(1/2)*r)/BesselI(0, (lambda/K_1)^(1/2)*a)

(18)

NULL

residual_at_pole1 := limit(integrand, lambda = 0)

1

(19)

# Pole at lambda = -K_1 alpha_n for n = 1,... N which are the N zeros of I_0(sqrt(a*lambda/K_1)) where alpha_n are the zeros of J_0(a*alpha_n)

 

NULL# One needs to assign the alpha_n to the roots of BesselI_0``


 

Download chaleur_1D_cylindrical_laplace.mw

@dharr 

I see you have now moved to Laplace on the other thread, which is how I would have done it

Yes I also think that this is the right way to go.

Here [1] means differentiate with respect to the first variable, and (a,t) is the point to avaluate at. For example 

Thank you, how would you specify second order derivative with respect to the first variable with this notation ? 

cheers,

Can

@Thomas Richard 

simplified_PDEL1 := eval(PDEL1,ic1);

Thank you, it was indeed a syntax error. Removing the brackets also worked.

And pdsolve should be used for PDEs; whereas dsolve is covering ODEs.

I am using Laplace transforms to end up with an ode to solve instead of the original pde.

Cheers,

@Rouben Rostamian  

Applying your suggestion to the above leads to

restart;
with(Physics);
vars := f(x);
sys_ode := diff(xhat(vars), f(x)) = a_1 + a_3*x, diff(yhat(vars), f(x)) = a_2 - a_4*y + f(x), diff(uhat(vars), f(x)) = a_3*u + (2*u)*a_4, diff(vhat(vars), f(x)) = a_4*v + u*diff(f(x), x);
                   
initvars := 0;
ics := xhat(initvars) = x, yhat(initvars) = y, uhat(initvars) = u, vhat(initvars) = v;
                      
dsolve({ics, eval(sys_ode)}, [xhat(f(x)), yhat(f(x)), uhat(f(x)), vhat(f(x))]);

Error, (in dsolve) found the following equations not depending on the unknowns of the input system: 
{uhat(0) = u, vhat(0) = v, xhat(0) = x, yhat(0) = y, 
(D(uhat))(f(x)) = a_3*u+2*a_4*u, 
(D(vhat))(f(x)) = a_4*v+u*(diff(f(x), x)), 
(D(xhat))(f(x)) = a_3*x+a_1, 
(D(yhat))(f(x)) = a_2-a_4*y+f(x)}

If I dont charge the module physics then I cannot go futher than

vars := f(x);
sys_ode := diff(xhat(vars), f(x)) = a_1 + a_3*x, diff(yhat(vars), f(x)) = a_2 - a_4*y + f(x), diff(uhat(vars), f(x)) = a_3*u + (2*u)*a_4, diff(vhat(vars), f(x)) = a_4*v + u*diff(f(x), x);

as I get

Error, invalid input: diff received f(x), which is not valid for its 2nd argument

Any idea on how to do ?

Cheers,

Can

@Rouben Rostamian  

Thank you for your help.

Is there a reason for looking at xhat as a function of all those variables?  Why not take xhat to be a function of a_1 only.  The solution will involve the remaining variables as parameters.

Yes, along the example given before, I was also trying to solve the following with maple:

with(Physics);
vars := x, y, u, v, a_1, a_2, a_3, a_4, f(x);
sys_ode := diff(xhat(vars), f(x)) = a_1 + a_3*x, diff(yhat(vars), f(x)) = a_2 - a_4*y + f(x), diff(uhat(vars), f(x)) = a_3*u + (2*u)*a_4, diff(vhat(vars), f(x)) = a_4*v + u*diff(f(x), x);
          vars := x, y, u, v, a_1, a_2, a_3, a_4, f(x)
initvars := x, y, u, v, a_1, a_2, a_3, a_4, 0;
ics := xhat(initvars) = x, yhat(initvars) = y, uhat(initvars) = u, vhat(initvars) = v;
dsolve({ics, eval(sys_ode)}, [xhat, yhat, uhat, vhat]);

it's an ode with respect to f(x) treaten as a derivative variable but f(x) and its diff(f,x) (its  derivative with respect to x) appears also on the RHS.

The derivative with respect to f(x) is writen as D_9 (as it differentiates with respect to the 9th argument of the variable-list) and does not match the argument f(x) that dsolve seems to be looking for. Hence (I guess) the error:

Error, (in dsolve) found the following equations not depending on the unknowns of the input system:
 {uhat(x, y, u, v, a_1, a_2, a_3, a_4, 0) = u, 
vhat(x, y, u, v, a_1, a_2, a_3, a_4, 0) = v, 
xhat(x, y, u, v, a_1, a_2, a_3, a_4, 0) = x, 
yhat(x, y, u, v, a_1, a_2, a_3, a_4, 0) = y, 
(D[9](uhat))(x, y, u, v, a_1, a_2, a_3, a_4, f(x)) = a_3*u+2*a_4*u, 
(D[9](vhat))(x, y, u, v, a_1, a_2, a_3, a_4, f(x)) = a_4*v+(diff(f(x), x))*u, 
(D[9](xhat))(x, y, u, v, a_1, a_2, a_3, a_4, f(x)) = a_3*x+a_1, 
(D[9](yhat))(x, y, u, v, a_1, a_2, a_3, a_4, f(x)) = a_2-a_4*y+f(x)}

 

Page 1 of 1