Maple 2022 Questions and Posts

These are Posts and Questions associated with the product, Maple 2022

with(Physics)

diff(x(t), `$`(t, 2)) = a(t)

diff(diff(x(t), t), t) = a(t)

(1)

dsolve(diff(diff(x(t), t), t) = a(t), arbitraryconstants = subscripted)

x(t) = Int(Int(a(t), t), t)+c__1*t+c__2

(2)

a(t) = 1.*Unit('m'/'s'^2), c__1 = 2*Unit('m'/'s'), c__2 = 3*Unit('m')

a(t) = 1.*Units:-Unit(m/s^2), c__1 = 2*Units:-Unit(m/s), c__2 = 3*Units:-Unit(m)

(3)

subs(a(t) = 1.*Units:-Unit(m/s^2), c__1 = 2*Units:-Unit(m/s), c__2 = 3*Units:-Unit(m), x(t) = Int(Int(a(t), t), t)+c__1*t+c__2)

x(t) = Int(Int(1.*Units:-Unit(m/s^2), t), t)+2*Units:-Unit(m/s)*t+3*Units:-Unit(m)

(4)

value(%)

x(t) = (1/2)*t^2*Units:-Unit(m/s^2)+2*Units:-Unit(m/s)*t+3*Units:-Unit(m)

(5)

NULL

Warning, units problem, not enough information to unambiguously deduce the units of the variables {t}; proceeding as if dimensionless

How can I provide the information to Maple that the unit of t is s?

 

Update: Temporarily disabling the warning would also be an option for the above case.

 

Download Unit_of_t.mw

I have 5 procedures saved in my library. Say Proc1, Proc2, Proc3, Proc4 and Proc5. They all work ok.

Proc5 calls Proc4 which calls Proc3 which calls Proc2 which calls Proc1. That is fine if I read them all in to the dodument.

Is it ok ok to put read commands in procedures 5..2 so I don't have to manually load them all myself into a new document or is there a better way to do this. These procedures are not part of a package.
I hace set 

currentdir();
      "C:\Users\Ronan\Documents\maple\Alibrary\Procedures"

which is where the procedures are stored. If I just load Proc5 it doesn't find the other procedures

dear all,

I am trying to use the solve to get symbolic answer.

to be simple, here is issue:

I have 15equations with 16 varibles.

the 15 equation as follows:

first 4 equation set 

equList_Ya :=

[ Va1 Ya11 + Va2 Ya12 + Va3 Ya13 + Va4 Ya14 = Ia1, 

  Va1 Ya21 + Va2 Ya22 + Va3 Ya23 + Va4 Ya24 = Ia2, 

  Va1 Ya31 + Va2 Ya32 + Va3 Ya33 + Va4 Ya34 = Ia3, 

  Va1 Ya41 + Va2 Ya42 + Va3 Ya43 + Va4 Ya44 = Ia4]

2nd 4 equation set

equList_Yb := 

[ Vb1 Yb11 + Vb2 Yb12 + Vb3 Yb13 + Vb4 Yb14 = Ib1, 

  Vb1 Yb21 + Vb2 Yb22 + Vb3 Yb23 + Vb4 Yb24 = Ib2, 

  Vb1 Yb31 + Vb2 Yb32 + Vb3 Yb33 + Vb4 Yb34 = Ib3, 

  Vb1 Yb41 + Vb2 Yb42 + Vb3 Yb43 + Vb4 Yb44 = Ib4]

the 3rd 4 equation set

equList_ConnectCondition :=

[Ia3 = -Ib1, Ia4 = -Ib2, Va3 = Vb1,  Va4 = Vb2]

the 4th 3 equation set

 equList_Y11Condition := [Va2 = 0, Vb3 = 0, Vb4 = 0]

so total equation number is 4+4+4+3=15

equList_Y11all := {Ia3 = -Ib1, Ia4 = -Ib2, Va2 = 0, Va3 = Vb1, Va4 = Vb2, Vb3 = 0, 
 Vb4 = 0, Va1*Ya11 + Va2*Ya12 + Va3*Ya13 + Va4*Ya14 = Ia1,
 Va1*Ya21 + Va2*Ya22 + Va3*Ya23 + Va4*Ya24 = Ia2,
 Va1*Ya31 + Va2*Ya32 + Va3*Ya33 + Va4*Ya34 = Ia3, 
 Va1*Ya41 + Va2*Ya42 + Va3*Ya43 + Va4*Ya44 = Ia4,
 Vb1*Yb11 + Vb2*Yb12 + Vb3*Yb13 + Vb4*Yb14 = Ib1,
 Vb1*Yb21 + Vb2*Yb22 + Vb3*Yb23 + Vb4*Yb24 = Ib2, 
 Vb1*Yb31 + Vb2*Yb32 + Vb3*Yb33 + Vb4*Yb34 = Ib3,
 Vb1*Yb41 + Vb2*Yb42 + Vb3*Yb43 + Vb4*Yb44 = Ib4}

the 16 varible is 

var_Y11 :=

[ Va1, Va2, Va3, Va4,

  Ia1, Ia2, Ia3, Ia4,

  Vb1, Vb2, Vb3, Vb4,

  Ib1, Ib2, Ib3, Ib4]

I hope to get Ia1=f(Va1) to calculate Y11(=Ia1/Va1)

but the solve give nothing:

Y11_all := solve(equList_Y11all, {Ia1}, symbolic = true);                        

Y11_all := NULL

Question1: why the above solve output nothing???

since not work,so I try another way:

polys_Y11 := map(lhs - rhs = 0, equList_Y11all);
    polys_Y11 := {Va2 = 0, Vb3 = 0, Vb4 = 0, Ia3 + Ib1 = 0, 

      Ia4 + Ib2 = 0, Va3 - Vb1 = 0, Va4 - Vb2 = 0, 

      Va1 Ya11 + Va2 Ya12 + Va3 Ya13 + Va4 Ya14 - Ia1 = 0, 

      Va1 Ya21 + Va2 Ya22 + Va3 Ya23 + Va4 Ya24 - Ia2 = 0, 

      Va1 Ya31 + Va2 Ya32 + Va3 Ya33 + Va4 Ya34 - Ia3 = 0, 

      Va1 Ya41 + Va2 Ya42 + Va3 Ya43 + Va4 Ya44 - Ia4 = 0, 

      Vb1 Yb11 + Vb2 Yb12 + Vb3 Yb13 + Vb4 Yb14 - Ib1 = 0, 

      Vb1 Yb21 + Vb2 Yb22 + Vb3 Yb23 + Vb4 Yb24 - Ib2 = 0, 

      Vb1 Yb31 + Vb2 Yb32 + Vb3 Yb33 + Vb4 Yb34 - Ib3 = 0, 

      Vb1 Yb41 + Vb2 Yb42 + Vb3 Yb43 + Vb4 Yb44 - Ib4 = 0}

and solve it with 

what I expected is Ia1=f(Va1) which Va1 as indepentent,but the above solution have no Va1 
Question2: is the (22) give the right answer for Ia1 ? what is the indepent var for the above solution?

here is my full maple worksheet

CoupleLineSeriesConnect_Question.mw

thx for your help.

best regards

Hello:

I've changed my laptop this week and I'd like to know how to install Maple with my account in the new one.

Thanks in advance.

Sergio Sanz.

I have some maple procedures that worked fine for many years. However, I started to notice some important performance deterioration with the new versions (Maple 2021 and 2022). I could identify the "culprit", they are some simplify/sum commands like this one:

A := simplify( sum( sum( sum( sum(sum( sum( sum( sum(
((x-i1)^2+(y-i2)^2)*((x-j1)^2+(y-j2)^2)*((x-k1)^2+(y-k2)^2)*((x-l1)^2+(y-l2)^2) ,i1=1..N)
,j1=1..N),k1=1..N),l1=1..N) ,i2=1..N),j2=1..N),k2=1..N),l2=1..N) );


Using, for instance, Maple 2017, this command requires to complete:

N:=4
memory used=258.2MB, alloc=44.3MB, time=1.58
N:=6
memory used=5247.1MB, alloc=394.0MB, time=39.35
N:=7
memory used=16405.4MB, alloc=775.9MB, time=135.78

On the other hand, in Maple 2022 (similar results for 2021), we have

N:=4
memory used=527.7MB, alloc=44.3MB, time=3.85
N := 6
memory used=10763.4MB, alloc=409.6MB, time=97.22
N := 7
memory used=34564.2MB, alloc=1139.9MB, time=351.94

My laptop (16G RAM) cannot execute N:=8 in Maple 2022 in reasonable time, I guess due to RAM limitation.

Does someone have any hint on how to improve the performance of these commands in Maple newer versions?
 

Mathmatica give the right answer:

I am composing a note using a maple document to demonstrate how the  mathematics from a text (that I find rather dense) can be unfolded into actual calculations using Maple.   My problem is to relate a typeset expression to actual Maple notation that bears no resemblance to the formal expression.

What I would like to do is to convert the typeset notation into a maple name which I can then associate with the calculation written in the functional notation.  For example, the quantity

|A0|2 := leftcontract(reverse(A0),A0);

When I use convert(|A0|2, name);   The result is abs(A[0] ^2 ) which isn’t typeset.

Similarly, the use of single left quotes around the expression as in ` |A0|2 ` it produces an error message when I test the result.  

type( ` |A0|2 ` ,name);
 Error, mismatched or missing bracket/operator

Admittedly, for ordinary calculations, just using any convenient name is more than satisfactory. However, in this case, where I’m trying to link together formal mathematical notation with the Maple operations it implies, being able to link the calculations to a special typeset same would be useful.

Let's say we have a procedure with an optional keyword option and a second procedure is calling this procedure and we want to give the user the option to set the kwarg of the first procedure in the second procedure as well. A simple example (just for the sake of the question, nothing meaningful in this example) is given below. Look at the kwarg "b" in test1. test2 is calling test1 and we want to have the option of setting "b" of test1 in test2 as well. But if I use "b = b" when calling test1, it doesn't work! I thought of using "`b` = b" and even "'b' = b", but they don't work either. One solution is to use a new name, say "c" and calling test1 by "b = c". But that is a bad choice. Because if you call test1 in so many other procedures, then you have to use so many names for one parameter, clearly this is not user friendly too, the user would prefer to remember a parameter by a fixed name. Is there any solution so that I can use the same name here?

test1 := proc( a :: posint, { b :: posint := 1 } ) :: posint:
	return( a + b ):
end proc:
test2 := proc( a :: posint, { b :: posint := 1 } ) :: posint:
	return( a * test1( a, b = b ) ):
end proc:
test3 := proc( a :: posint, { b :: posint := 1 } ) :: posint:
	return( a * test1( a, `b` = b ) ):
end proc:
test4 := proc( a :: posint, { c :: posint := 1 } ) :: posint:
	return( a * test1( a, b = c ) ):
end proc:

I use the command line mint  since all my code in .mpl files.

First question:

I noticed mint complains that module name is global, for proc inside the module itself, when the name is used as type of a local variable. But maplemint does not complain. Here is an example

A:=module()
   local B:=module()
         option object;
         export n::integer:=1;
   end module;

   export foo:=proc()
      local a::A:-B;  #mint complains that A is global not declared!
      a:=Object(A:-B);
      a:-n:=2;
   end proc;
end module;

running mint -i 2 A.mpl gives

These names were used as global names but were not declared:  A

But A is the module name where the whole code is sitting inside it?  To fix this, I have to add

A:=module()
   local B:=module()
         option object;
         export n::integer:=1;
   end module;

   export foo:=proc()
      global A;  ###################add this
      local a::A:-B;
      a:=Object(A:-B);
      a:-n:=2;
   end proc;
end module;

And now the message/warning is gone. But the above looks really strange.  maplemint does not complain about global A:

But notice that maplemint gives waring that n is never used but mint do not. Another difference!

Which is right about the global name message?

Second question is: I do not understand the message

     These local variables were assigned a value, but otherwise unused:  B

which both mint and maplemint give,

What does it mean B is not used?? I used it to make object a inside proc foo().

What would one do to get rid of this message?

Final comment:

I find many message come out that are not real problems at all. For example, if I declare local variable and use it in equation, as symbol, mint complains that the variable was not assigned a value. Here is an example

foo:=proc()
   local x,eq;
   eq := x^2;
   return eq;
end proc;

now maplemint(foo) gives

Procedure foo()
  These local variables were used but never assigned a value:  x

I know there are ways to filter out these messages. But I am afraid if I do that, I will filter out a message that indicates a real problem?

According to help   -i 2 shows Display severe and serious errors (default)

if I use -i 1, then these message do not show. What do others use?  Level 1 or 2? If I use level 3, then more strange messages show which I do not understand at all how to remove. Such as

These parameters or local variables are also system defined names:
      thismodule
  These names are special to Maple:
      thismodule

So I stick to level 2 for now. but 90% of the messages I get are not real errors at all.

I'm trying to get the plot command to display correctly on a logarithmic axis. The calculations in Maple are correct but the plots using logarithmic axis are wrong and start at different points on the x-axis. Is there a way to force Maple to display set points along the axis? e.g. 1, 3, 5, 10 of each logarithmic range. 

Log_plot_example.mw

When running this

expr:=5;
for item in expr do
    print(item);
od;

it prints 5 as expected. When when running

expr:=I;  (* NOTICE this is complex  I not the number 1 *)
for item in expr do
    print(item);
od;

it prints the number 1 and not I

I am sure there is a good reason why this happens.

But what should one do to insure they get the complex I in the second example when iterating over a sum of numbers, which in this example happened to be just one number who is complex? I know I can add explicit check to avoid this edge case.

The problem is that I want to iterate of sum of iterms, One or more of them can be complex I. (it is a result of doing series expansion of function at infinity, and I want to iterate over each term in the series one by one).

expr:=9+I:
for item in expr do
    print(item);
od;

Gives   9,1   and not 9,I

While

expr:=a+I:
for item in expr do
    print(item);
od;

does now give exected output   I,a  and not 1,a

 

What should one do to insure the for loop always get each term in the sum, even if it is complex?

I can't check the item inside the loop, because by then it is too late as the compelx I is lost already.

Is there a better way to iterate over sum of terms and look at each term as is and not lose the complex I in the way?

Maple 2022.1

Update 

For now, I am doing this hack. Since the input is always a sum (with + or - terms), then I convert the input to string, split on delimitors and then parse the entries back to Maple and now it is a list so I can iterate over each term. 

expr:=I+3:
String(expr):
StringTools:-SubstituteAll(%,"+","|+"):
StringTools:-SubstituteAll(%,"-","|-"):
StringTools:-Split(%,"|"):
map(Z->parse(Z),%):
for item in % do
    print(item);
od;

And when expr is

expr:=1+5*I+sin(x)-sqrt(8)*I-a;

It is a hack, but I could not find a relaible way to iterate over each term in the sum of terms without losing the complex I if one term was complex. I need to test this more. It is meant to work only on expression which is sum of terms, which is where I will use it.

Compare the output of the above for expr:=1-I;

The for loop gives 1,-1  but the string hack gives 1,-I which is what I wanted.

Ofcourse the above could fail if there is a "-" or "+" inside the term itself (for example  1+I+sin(1-x) , but this do not happen for the cases I am using this for, which just looking at terms of series expansion around either zero or infinity. These will be just power series terms separated by "+" or possibly "-"

If there is better way do this, that will be great.

 

 

 

Given a module A, it has a proc which is called to set some internal variable to some value. This proc is meant to be used as initialization of the module, and not meant to return anything (i.e. procedure vs. a function in other languages, where a procedure does not return anything, but a function does).

When this proc is called from outside, the result of the last assignment in the proc is returned back to the user since that is the default behavior.

Is there a way to prevent this, other than adding an explicit NULL at the end of this proc?  Adding : at the end the last statement or at the end of the proc did not prevent this.

Here is an example

restart;
A:=module()
  local r::integer;
  export set_r:=proc(r::integer)   
    A:-r:= r:   
  end proc:
end module:

Now when calling  A:-set_r(10), Maple will return back/echo back 10. But this is not something I want.

A:-set_r(10);

    10

To prevent this, currently I add NULL; at end of the proc, like this

A:=module()
  local r::integer;
  export set_r:=proc(r::integer)   
    A:-r:= r:   
    NULL;
  end proc:
end module:

and now

A:-set_r(10);

returns nothing. Well, it returns NULL but that is like nothing.

My question is, is there a better way to do this? Can one define a proc in Maple that returns nothing?

Maple 2022.1

I want to replace   expression that contains function whose first argument is always R0 with r0.

I can do it when the function is  f(R0) or f(R0,x) but I do not know how to tell Maple to do the replacement of the first argument, if the function has 0 or more arguments beyond R0.  This is what I tried

restart;

#this works
expr:=A+f(R0);
evalindets(expr, 'anyfunc(identical(R0))', Z-> subsop(1 = r0, Z ));

A+f(R0)

A+f(r0)

#this also works, when there is one argument after R0
expr:=A+f(R0,x);
evalindets(expr, 'anyfunc(identical(R0),anything)', Z-> subsop(1 = r0, Z ));

A+f(R0, x)

A+f(r0, x)

#but how to do it for any number of arguments after R0? This does not work
expr:=A+f(R0,x,y);
evalindets(expr, 'anyfunc(identical(R0),anything)', Z-> subsop(1 = r0, Z ));

A+f(R0, x, y)

A+f(R0, x, y)

#I do not want to keep adding anything to match the number of argument, as I do not know how many there could be
expr:=A+f(R0,x,y);
evalindets(expr, 'anyfunc(identical(R0),anything,anything)', Z-> subsop(1 = r0, Z ));

A+f(R0, x, y)

A+f(r0, x, y)

 

Download oct_11_2022_anyfunc.mw

I do not want to keep writing anything,anything,anything, for as many times as there could be more arguments.  I wan to say something like

evalindets(expr, 'anyfunc(identical(R0),anything_zero_or_more_arguments)', Z-> subsop(1 = r0, Z ));

Is there a better way to do this whole thing than what I am doing?

Update

Just after I posted this question, I found another structured type that worked

expr:=A+f(R0,x,y,z,h);
evalindets(expr, 'patfunc(identical(R0),anything)', Z-> subsop(1 = r0, Z ));

So patfunc works for zero or more argument after R0, which is what I wanted.

Given that types are core part of Maple, I find it so amazing that the help page "Definition of a Structured Type in Maple"  has so little examples and description of all these types in it. I can't figure what 80% of the types there do, since there are no examples.

Maplesoft should really do a better job on this page and add much much more examples and description. After all, strong typing is what is supposed to be the strength of Maple. But looking at this help page, it does not give one this impression at all.

I have a matrix and generate difference tables from it's diagonals. Works fine. Then I make a polynomial from the first row of the difference table. I have three polynomials. They are all acting inert, just will not evaluate.

I cannot see what the problem is.
 

restart

interface(rtablesize = 12)

[10, 10]

i) Recurrence relation

 

"Recurrence Relation for next column. . The matrix starts at 1,1 not 0,0 A(i,j)=A(i-1,j-1)*(j-2)+A(i,j-1)"
``

NULL``

Msize := 12

A := Matrix(Msize)

A[1, 1] := 1

1

for i from 2 to Msize do for j from i to Msize do A[i, j] := A[i, j-1]*(j-2)+A[i-1, j-1] end do end do

A

Matrix(%id = 36893489695842121836)

``

``

viia)  General formula for diagonals:-

Select diagonal  & Difference Table procedures

 

diag := proc (M::Matrix, dg::posint) local i, rd, c, dt; rd := LinearAlgebra:-RowDimension(M); dt := Matrix(rd, rd); for i from dg to rd do dt[1+i-dg, 1] := M[1+i-dg, i] end do; return dt end proc

proc (M::Matrix, dg::posint) local i, rd, c, dt; rd := LinearAlgebra:-RowDimension(M); dt := Matrix(rd, rd); for i from dg to rd do dt[1+i-dg, 1] := M[1+i-dg, i] end do; return dt end proc

diag(A, 5)

Matrix(%id = 36893489695866397388)

diftab := proc (M::Matrix) local i, j, cl, rd; rd := LinearAlgebra:-RowDimension(M); for i from rd by -1 to 1 do if M[i, 1] = 0 then rd := rd-1 else break end if end do; cl := rd; for i from 2 to cl do for j to rd-i do M[j, i] := M[j+1, i-1]-M[j, i-1] end do end do; print(M); return M end proc

proc (M::Matrix) local i, j, cl, rd; rd := LinearAlgebra:-RowDimension(M); for i from rd by -1 to 1 do if M[i, 1] = 0 then rd := rd-1 else break end if end do; cl := rd; for i from 2 to cl do for j to rd-i do M[j, i] := M[j+1, i-1]-M[j, i-1] end do end do; print(M); return M end proc

diftab(diag(A, 4))

Matrix(%id = 36893489695866348116)

gnrpoly := proc (M::Matrix) local i, p, rd, n; rd := LinearAlgebra:-RowDimension(M); for i from rd by -1 to 1 do if M[1, i] = 0 then rd := rd-1 else break end if end do; p := 0; for i from 0 to rd-1 do p := p+M[1, i+1]*binomial(n, i) end do; return p end proc

proc (M::Matrix) local i, p, rd, n; rd := LinearAlgebra:-RowDimension(M); for i from rd by -1 to 1 do if M[1, i] = 0 then rd := rd-1 else break end if end do; p := 0; for i from 0 to rd-1 do p := p+M[1, i+1]*binomial(n, i) end do; return p end proc

"for i from 0 to 3 do print("Diagonal  ",i);   eqn:=gnrpoly(diftab(diag(A,i+1)));    cat(poly,i):=unapply(factor(expand(eqn)),n);    print("--------------------------------------------------------------");  end do"

"--------------------------------------------------------------"

NULL

``

``

The equations do not evaluate. They are acting inert

poly1(n)

(1/2)*n*(1+n)

poly1(2)

(1/2)*n*(1+n)

poly2(n)

(1/24)*n*(3*n+5)*(n+2)*(n+1)

poly2(3)

(1/24)*n*(3*n+5)*(n+2)*(n+1)

poly3(n)

(1/48)*n*(n+1)*(n+3)^2*(n+2)^2

poly3(-1)

(1/48)*n*(n+1)*(n+3)^2*(n+2)^2

eqn

6*n+38*binomial(n, 2)+93*binomial(n, 3)+111*binomial(n, 4)+65*binomial(n, 5)+15*binomial(n, 6)

value(eval(eqn, n = 4))

6*n+38*binomial(n, 2)+93*binomial(n, 3)+111*binomial(n, 4)+65*binomial(n, 5)+15*binomial(n, 6)

eqn1 := expand(eqn)

(3/4)*n+2*n^2+(97/48)*n^3+(47/48)*n^4+(11/48)*n^5+(1/48)*n^6

eval(eqn1, [n = 4])

(3/4)*n+2*n^2+(97/48)*n^3+(47/48)*n^4+(11/48)*n^5+(1/48)*n^6

NULL


 

Download Q_10-10-22_gentrate_eqn.mw

Hi,

I have developed a series of two dimension PDE equations (x,t) and I am going to solve them numerically in Maple. 

However, after putting all equations, initial conditions and boundary conditions, I faced couple of errors. I would be grateful if you kindly look at the code and let me know how I can solve it.

Thank you

Here is the code:


restart: with(plots): with(LinearAlgebra):

L:=0.003:   # thickness
rho_w:=997.77:
rho_s:=1419:
lambda_s:=0.46:
lambda_w:=0.644:
lambda_g:=0.026:
cp_s:=3734:
cp_w:=4183:
cp_v:=1900:
cp_a:=1005.68:
M_v:=18.016:
M_a:=28.966:
R:=8.314:
epsilon:= t -> 0.9*(1-t/10):
Cw0:=6:
Sw0:= Cw0/rho_w/epsilon(0):
pi:=3.1415:
p0:=patm:
patm:=10^5:
T0:=256:
p_air:=0.2*19724:
h_m:=0.017904:
h_T:=16.746:
T_air:=380:

Xdb:= (x,t) -> S(x,t):
Cw:= (x,t) ->  rho_w*epsilon(t)*S(x,t):  
Cg:= (x,t) -> rho_g(x,t)*epsilon(t)*(1-S(x,t)):
Cv:= (x,t) -> rho_v(x,t)*epsilon(t)*(1-S(x,t)):
Ca:= (x,t) -> rho_a(x,t)*epsilon(t)*(1-S(x,t)):
nw:= (x,t) -> -rho_w* k_rw(x,t)*K(t)/(mu_w(x,t))*(diff(p(x,t),x)-D_c(x,t)*diff(Cw(x,t),x) ):
ng:= (x,t) -> -rho_g(x,t)* k_rg(x,t)*K(t)/(mu_g(x,t))*(diff(p(x,t),x)):
nv:= (x,t) -> -w_v(x,t)*rho_g(x,t)* k_rg(x,t)*K(t)/(mu_g(x,t))*(diff(p(x,t),x))-binary(x,t):
na:= (x,t) -> -(1-w_v(x,t))*rho_g(x,t)* k_rg(x,t)*K(t)/(mu_g(x,t))*(diff(p(x,t),x))+binary(x,t):
M_g:= (x,t) -> M_v*w_v(x,t)+M_a*(1-w_v(x,t)):
rho_g:= (x,t) -> p(x,t)*M_g(x,t)/R/T(x,t):
rho_v:= (x,t) -> rho_g(x,t)*w_v(x,t):
rho_a:= (x,t) -> rho_g(x,t)*(1-w_v(x,t)):
binary:= (x,t) -> rho_g(x,t)*epsilon(t)*(1-S(x,t))*Deff(x,t)*diff(w_v(x,t),x):
Deff:= (x,t) -> 2.3*10^(-5)*p0/p(x,t)*(T(x,t)/T0)^1.81*(epsilon(t)*(1-S(x,t)))^(4/3):
mu_w:= (x,t) -> rho_w*exp(-19.143+1540/T(x,t)):
mu_g:= (x,t) -> 0.017/1000*(T(x,t)/273)^0.65:
p_veq:= (x,t) -> p_vsat(x,t)*a_w(x,t):
p_vsat:= (x,t) -> exp(-5800.2206/T(x,t)+1.3915-0.0486*T(x,t)+0.4176*10^(-4)*T^2-0.01445*10^(-7)*T^3+6.656*ln(T(x,t))):
a_w:= (x,t) -> exp(-0.182*Xdb(x,t)^(-0.696)+0.232*exp(-43.949*Xdb(x,t))*Xdb(x,t)^0.0411*ln(p_vsat(x,t))):
h_fg:= (x,t) -> 3167.2-2.432*T(x,t):
I_vap:= (x,t) -> M_v*K_eff*(p_veq(x,t)-p(x,t))/R/T(x,t):
K_eff:=1000:
rhocp:= (x,t) -> w_v(x,t)*epsilon(t)*(1-S(x,t))*rho_g(x,t)*cp_v+(1-w_v(x,t))*epsilon(t)*(1-S(x,t))*rho_g(x,t)*cp_a+epsilon(t)*S(x,t)*rho_w*cp_w+(1-epsilon(t))*rho_s*cp_s:
lambda:= (x,t) -> epsilon(t)*(1-S(x,t))*lambda_g+epsilon(t)*S(x,t)*lambda_w+(1-epsilon(t))*lambda_s:
ncp:= (x,t) -> nv(x,t)*cp_v+na(x,t)*cp_a+nw(x,t)*cp_w:
k_rw:= (x,t) -> S(x,t)^3:
k_rg:= (x,t) -> 1.01*exp(-10.86*S(x,t) ):
K:= t -> 10^(-10)*(1-t/10):
D_c:= (x,t) -> 10^(-8)*exp(-2.8+2*Xdb(x,t)):

PDE_m1:=diff(Cw(x,t),t)+ diff(nw(x,t),x)=-I_vap(x,t) :
PDE_m2:=diff(Ca(x,t),t)+ diff(na(x,t),x)=0 :
PDE_m3:=diff(Cg(x,t),t)+ diff(ng(x,t),x)=I_vap(x,t) :
PDE_T:=diff(rhocp(x,t),t)+ diff(ncp(x,t)*T(x,t),x)=diff( lambda(x,t)*diff(T(x,t),x),x) -h_fg(x,t)*I_vap(x,t):

IBC_S:={S(x,0)=Sw0, D[1](Cw)(0,t)=0, subs(x=L,nw(x,t))=epsilon(t)*S(L,t)*h_m/R/T(L,t)*(p_veq(L,t)-p_air) }:
IBC_p:={D[1](p)(0,t)=0,p(L,t)=patm,p(x,0)=patm}:
IBC_T:={h_T*(T(L,t)-T_air)+epsilon(t)*S(L,t)*h_m/R/T(L,t)*(p_veq(L,t)-p_air)*h_fg(L,t)=0, T(x,0)=T_air, D[1](T)(0,t)=0}:
IBC_w:={w_v(x,0)=0.0262, subs(x=L,nv(x,t))=epsilon(t)*(1-S(L,t))*h_m/R/T(L,t)*(p_veq(L,t)-p_air), D[1](w_v)(0,t)=0 }:

pds:=pdsolve(PDE_m1,PDE_m2,PDE_m3,PDE_T,{IBC_S,IBC_p,IBC_T,IBC_w},numeric);
 

First 30 31 32 33 34 35 36 Last Page 32 of 43