mmcdara

7891 Reputation

22 Badges

9 years, 58 days

MaplePrimes Activity


These are replies submitted by mmcdara

@mehran rajabi 

Could be interested by this?
 


 

restart;

DiffGivenVars := proc(f, given)
   local vars, u:
   vars := indets(f):
   u    := vars minus given:
   diff(f, op(u))[op(given)]
end proc:

Warning, `u` is implicitly declared local to procedure `DiffGivenVars`

 

f := (P, V, T) ->  P*V/T

proc (P, V, T) options operator, arrow; P*V/T end proc

(1)

DiffGivenVars(f(P, V, T), {P, T});
DiffGivenVars(f(P, V, T), {P});

(P/T)[P, T]

 

(-P/T^2)[P]

(2)

# In case your equation of states would depend on some formal constants

DiffGivenVars2 := proc(f, vars, WithRespectTo)
   local u:
   u := vars minus WithRespectTo:
   diff(f, op(WithRespectTo))[op(u)]
end proc:

f := V^3-R*T*V^2/P-(B^2+P*B*R*T*sqrt(T)/(P-A))*V-A*B/(P*sqrt(T))

V^3-R*T*V^2/P-(B^2+P*B*R*T^(3/2)/(P-A))*V-A*B/(P*T^(1/2))

(3)

DiffGivenVars2(f, {P, T, V}, {V});
DiffGivenVars2(f, {P, T, V}, {T, P})

(3*V^2-2*R*T*V/P-B^2-P*B*R*T^(3/2)/(P-A))[P, T]

 

(R*V^2/P^2-((3/2)*B*R*T^(1/2)/(P-A)-(3/2)*P*B*R*T^(1/2)/(P-A)^2)*V-(1/2)*A*B/(P^2*T^(3/2)))[V]

(4)

g := V^3-R*T*V^2/P-(B^2+P*B*R*T*sqrt(T)/(P-A))*V-A*B/(P*sqrt(T)) = 0:

pressure := [ solve(g, P) ]

[(1/2)*(-T^(3/2)*R*V^2+A*B^2*T^(1/2)*V-A*T^(1/2)*V^3-B*A+(4*A*B*T^(7/2)*R^2*V^3+T^3*R^2*V^4+2*A*B^2*T^2*R*V^3-2*A*T^2*R*V^5+A^2*B^4*T*V^2-2*A^2*B^2*T*V^4+A^2*T*V^6+4*A^2*B^2*R*T^2*V+2*A*B*T^(3/2)*R*V^2+2*A^2*B^3*T^(1/2)*V-2*A^2*B*T^(1/2)*V^3+A^2*B^2)^(1/2))/(V*(B*R*T^2+B^2*T^(1/2)-T^(1/2)*V^2)), (1/2)*(-T^(3/2)*R*V^2+A*B^2*T^(1/2)*V-A*T^(1/2)*V^3-B*A-(4*A*B*T^(7/2)*R^2*V^3+T^3*R^2*V^4+2*A*B^2*T^2*R*V^3-2*A*T^2*R*V^5+A^2*B^4*T*V^2-2*A^2*B^2*T*V^4+A^2*T*V^6+4*A^2*B^2*R*T^2*V+2*A*B*T^(3/2)*R*V^2+2*A^2*B^3*T^(1/2)*V-2*A^2*B*T^(1/2)*V^3+A^2*B^2)^(1/2))/(V*(B*R*T^2+B^2*T^(1/2)-T^(1/2)*V^2))]

(5)

# derivative of the pressure given the temperature

DiffGivenVars2(pressure[1], {T, V}, {V});

((1/2)*(-2*T^(3/2)*R*V+A*B^2*T^(1/2)-3*A*T^(1/2)*V^2+(1/2)*(12*A*B*T^(7/2)*R^2*V^2+4*T^3*R^2*V^3+6*A*B^2*T^2*R*V^2-10*A*T^2*R*V^4+2*A^2*B^4*T*V-8*A^2*B^2*T*V^3+6*A^2*T*V^5+4*A^2*B^2*R*T^2+4*A*B*T^(3/2)*R*V+2*A^2*B^3*T^(1/2)-6*A^2*B*T^(1/2)*V^2)/(4*A*B*T^(7/2)*R^2*V^3+T^3*R^2*V^4+2*A*B^2*T^2*R*V^3-2*A*T^2*R*V^5+A^2*B^4*T*V^2-2*A^2*B^2*T*V^4+A^2*T*V^6+4*A^2*B^2*R*T^2*V+2*A*B*T^(3/2)*R*V^2+2*A^2*B^3*T^(1/2)*V-2*A^2*B*T^(1/2)*V^3+A^2*B^2)^(1/2))/(V*(B*R*T^2+B^2*T^(1/2)-T^(1/2)*V^2))-(1/2)*(-T^(3/2)*R*V^2+A*B^2*T^(1/2)*V-A*T^(1/2)*V^3-B*A+(4*A*B*T^(7/2)*R^2*V^3+T^3*R^2*V^4+2*A*B^2*T^2*R*V^3-2*A*T^2*R*V^5+A^2*B^4*T*V^2-2*A^2*B^2*T*V^4+A^2*T*V^6+4*A^2*B^2*R*T^2*V+2*A*B*T^(3/2)*R*V^2+2*A^2*B^3*T^(1/2)*V-2*A^2*B*T^(1/2)*V^3+A^2*B^2)^(1/2))/(V^2*(B*R*T^2+B^2*T^(1/2)-T^(1/2)*V^2))+(-T^(3/2)*R*V^2+A*B^2*T^(1/2)*V-A*T^(1/2)*V^3-B*A+(4*A*B*T^(7/2)*R^2*V^3+T^3*R^2*V^4+2*A*B^2*T^2*R*V^3-2*A*T^2*R*V^5+A^2*B^4*T*V^2-2*A^2*B^2*T*V^4+A^2*T*V^6+4*A^2*B^2*R*T^2*V+2*A*B*T^(3/2)*R*V^2+2*A^2*B^3*T^(1/2)*V-2*A^2*B*T^(1/2)*V^3+A^2*B^2)^(1/2))*T^(1/2)/(B*R*T^2+B^2*T^(1/2)-T^(1/2)*V^2)^2)[T]

(6)

 


 

Download derivative.mw

The name Pr cannot be at the same time on the left and on the right of the '=' sign.
I think you should clarify your question.

In the meantime tere are a fex possibilities:


 

restart:

Option 1

Pr := LinearAlgebra:-DiagonalMatrix(<-1, -2/3, -1/3, 0>);

Pr := Matrix(4, 4, {(1, 1) = -1, (2, 2) = -2/3, (3, 3) = -1/3, (4, 4) = 0}, storage = diagonal, shape = [diagonal])

(1)

# or more generaly :

Pr := x -> LinearAlgebra:-DiagonalMatrix(x):

x := <-1, -2/3, -1/3, 0, 8, -2>:
Pr(x)

Matrix([[-1, 0, 0, 0, 0, 0], [0, -2/3, 0, 0, 0, 0], [0, 0, -1/3, 0, 0, 0], [0, 0, 0, 0, 0, 0], [0, 0, 0, 0, 8, 0], [0, 0, 0, 0, 0, -2]])

(2)

Option 2

# the name "Pr" cannot be used in the definition of matrix Pr (recursivity conflict)

Pr := N -> LinearAlgebra:-DiagonalMatrix(<seq(f(x__||n), n=0..N-1)>):

Pr(6);

Matrix([[f(x__0), 0, 0, 0, 0, 0], [0, f(x__1), 0, 0, 0, 0], [0, 0, f(x__2), 0, 0, 0], [0, 0, 0, f(x__3), 0, 0], [0, 0, 0, 0, f(x__4), 0], [0, 0, 0, 0, 0, f(x__5)]])

(3)

# or, for cosmetic display only ...

Pr := N -> LinearAlgebra:-DiagonalMatrix(<seq('Pr'(x__||n), n=0..N-1)>):

'Pr' = Pr(6);

 

Pr = (Matrix(6, 6, {(1, 1) = Pr(`#msub(mi("x"),mi("0"))`), (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (1, 5) = 0, (1, 6) = 0, (2, 1) = 0, (2, 2) = Pr(`#msub(mi("x"),mi("1"))`), (2, 3) = 0, (2, 4) = 0, (2, 5) = 0, (2, 6) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = Pr(`#msub(mi("x"),mi("2"))`), (3, 4) = 0, (3, 5) = 0, (3, 6) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = Pr(`#msub(mi("x"),mi("3"))`), (4, 5) = 0, (4, 6) = 0, (5, 1) = 0, (5, 2) = 0, (5, 3) = 0, (5, 4) = 0, (5, 5) = Pr(`#msub(mi("x"),mi("4"))`), (5, 6) = 0, (6, 1) = 0, (6, 2) = 0, (6, 3) = 0, (6, 4) = 0, (6, 5) = 0, (6, 6) = Pr(`#msub(mi("x"),mi("5"))`)}))

(4)

Option 3 : but maybe Pr is a matrix made of bloc matrices ???

with(Student[LinearAlgebra]):

BlockSize := 2:
N := 3:

x := [seq(x__||n, n=0..N-1)];
Matrix(BlockSize$2, symbol=u):
blocks := [ seq(Matrix(BlockSize$2, (i,j) -> f[i,j](x[n])), n=1..N) ]:
DiagonalMatrix(blocks);

x := [`#msub(mi("x"),mi("0"))`, `#msub(mi("x"),mi("1"))`, `#msub(mi("x"),mi("2"))`]

 

Matrix([[f[1, 1](x__0), f[1, 2](x__0), 0, 0, 0, 0], [f[2, 1](x__0), f[2, 2](x__0), 0, 0, 0, 0], [0, 0, f[1, 1](x__1), f[1, 2](x__1), 0, 0], [0, 0, f[2, 1](x__1), f[2, 2](x__1), 0, 0], [0, 0, 0, 0, f[1, 1](x__2), f[1, 2](x__2)], [0, 0, 0, 0, f[2, 1](x__2), f[2, 2](x__2)]])

(5)

 


 

Download Pr.mw

@Carl Love 

 

For info, these "when V is constant" or "when T is constant" expressions, which often seem strange to "pure" mathematicians, are of common practice in thermophysics. 
Here is a reference, maybe not the best one, only one of the first I found on the web, which speaks about that
You will be able to find many others if you have a few time to spend on this topic (just type for instance "derivative of pressure with respect to volume" in your explorator)

https://www.harding.edu/lmurray/themo_files/notes/ch04.pdf

 

@vv 

INo, the question is general and noy limited to uniform distributions.

I tried to manipulate the inequality and I have found some interesting results (unless I did some mistake, a thing I'm not unfamiliar with).
I think I proved there exist, for any triple (a, b, c) a value m of n such that 
(a^(m+1) - b^(m+1)) - c * (a^m- b^m) > 0
and that, for any positive integer p, (a^(m+p+1) - b^(m+p+1)) - c * (a^m+p- b^m+p) > 0.

Here is the "pseudo worksheet", could you be kind enough to check it?
Maybe this could give some hint to provide elements for a "non statistical" answer.

Thanks in advance


 

restart:

# (a^(n+1) - b^(n+1)) - c * (a^n - b^n) =  a^n * (a-c) -  b^n * (b-c)

f := (a, b, c, n) -> a^n * (a-c) -  b^n * (b-c)

proc (a, b, c, n) options operator, arrow; a^n*(a-c)-b^n*(b-c) end proc

(1)

f(a, b, c, 0);  # always negative for b > a

a-b

(2)

# when n tends to + infinity the behavior of f is the one of
# - b^n * (b-c)
# then the limit of f (n --> + infinity) is positive
#
# prop1 : there exist an integer m > 0 such that f(a, b, c, m) > 0

# try to find the signd of f(a, b, c, m):

f1 := f(a, b, c, m+1);

# f1 := a * ( a^m*(a-c)-b^m*(b-c) ) + a*b^m*(b-c) - b*b^m*(b-c)
# f1 := a * ( a^m*(a-c)-b^m*(b-c) ) + b^m*(b-c)*(a-b)

# verify
expand(a * ( a^m*(a-c)-b^m*(b-c) ) + b^m*(b-c)*(a-b) -f1);



 

a^(m+1)*(a-c)-b^(m+1)*(b-c)

 

0

(3)

# by definition m is such that a^m*(a-c)-b^m*(b-c) > 0
#
# then a * ( a^m*(a-c)-b^m*(b-c) ) > 0
#
# on the other hand (b-c)*(a-b) > O because a < b < c
#
# thus b^m*(b-c)*(a-b)  and finally
#
# f(a, b, c, m) > 0 ==> f(a, b, c, m+1) > 0

# prop2 : Let m such that f(.., m) > 0 then, for any strictly positive integer p,
#         f(..., m+p) > 0

 


 

Download abc.mw

@vv @Preben Alsholm

Questionable and too simplistic a way to proceed.

Put "a" aside to simplify the problem and do visualize what happens in 2D.

What you did is
  1/ generate c in [0,1]
   2/ generate b in [0, c].
But you  could also have done this by switching steps 1 and 2:
  1/ generate b in [0,1]
   2/ generate c in [1, b].

Intuitively one could expect the results to be the same. But sampling B conditionnally to C (first way) is generally not the same thing than to samplong C condiitonnally to B (second way)

Please look to the attached file to see the differences between the two plots

Which means the result you give (7393 "OK" out of 10^5) depends on the way you sample the three uniform random variables  A, B and C whose outcomes are a, b anc c (last part of the attached file)

A famous example of such a problem  is given by Bertrand's paradox

 


First way to generate 0 < b < 1:

  1/ generate c : 0 < c < 1
  2/generate b : 0 < b < c

restart:

interface(rtablesize=10):

N := 10^4:
r := LinearAlgebra:-RandomMatrix(N, 2, generator=rand(0. .. 1.)):

bc := < r[..,1] *~ r[..,2] | r[.., 2] >:
# second column of bc represents c in [0,1]
# first column of bc now represent b in [0,c]

Statistics:-ScatterPlot(bc[..,1] , bc[..,2], labels=[b, c], symbol=point )

 


second way to generate 0 < b < 1:


  1/ generate b :  0 < b < 1
  2/ generate c : b < c < 1

bc := < r[.., 1] | r[.., 1] + (1 -~ r[.., 1]) *~ r[.., 2] >:
Statistics:-ScatterPlot(bc[..,1] , bc[..,2], labels=[b, c], symbol=point )

 

N := 10^5:
r := LinearAlgebra:-RandomMatrix(N, 3, generator=rand(0. .. 1.)):


c := r[.., 3]:
b := r[.., 2] *~ c:
a := r[.., 1] *~ b:

n  := 6: # for example

OK := N - add(Heaviside~((a^~(n+1) - b^~(n+1)) - c *~ (a^~n - b^~n) ));

HFloat(3100.0)

(1)

a := r[.., 1]:
b := a + (1 -~ a) *~ r[.., 2]:
c := b + (1 -~ b) *~ r[.., 3]:

OK := N - add(Heaviside~((a^~(n+1) - b^~(n+1)) - c *~ (a^~n - b^~n) ));

HFloat(35349.0)

(2)

 


 

Download BandC.mw

@Carl Love @acer

Than to both of you for your detailed answers 

@vv 

When those problems are challenges in mathematical olympiads, which could be the case of this one, the solution that is expected is more related to the "aha effect" mentionned by acer than to notions of advanced algebra.
That's why I sent this socratic.org link to LemonySnicket

But everybody is free to proposed any other standard method to solve the problem, and we have had a lot here

Kind of headache one can find in some mathematical problem books.

I guess you were interested in some astute method to find the solutions by hand?
If so look here 

https://socratic.org/questions/if-a-b-c-1-a-2-b-2-c-2-2-a-3-b-3-c-3-3-then-find-the-value-of-a-4-b-4-c-4

@dharr 

Thanks dharr, you make me feel a wit bit less dumb.

@Rouben Rostamian

Don't feel obliged to apologize because I didn't feel offended. As English is not my mother tongue, I also express myself in a way that may seem harsh.
I reluctantly analyze my code, sure I was it was correct.
Unfortunately I was wrong and it contained an error. In fact I did not impose the boundary conditions at all (funny thing given it was our main topic of debate)

Here is the revised version.
The guilty lines are red written on yellow background.
I will continue to review my code to ensure it contains no further errors. I'll inform you if it's not the case.

May we consider this exchange is over?


 

restart;

with(plots):

interface(rtablesize=10):

f     := unapply(-x^2+1,x):
mu[1] := t -> 0:
mu[2] := t -> 0:
g     := (t, x) -> 0:

l   := 2:
n   := 50:
h   := l/n:
T   := 0.5:
eps := 0.01:
CFL := 1/2-eps:
tau := CFL*h^2:
m   := round(T/tau):

CFL := tau/h^2;
if CFL >= 1/2 then WARNING("The CFL condition is not verified, the scheme is unstable") end if;

.4900000000

(1)

FD_T  := k -> (yF[k]-y[k]) / tau:
FD_X  := k -> (y[k+1]-y[k]) / h:
CDD_X := k -> simplify((FD_X(k,j)-FD_X(k-1,j)) / h):

eq := (k,j) -> FD_T(k) = CDD_X(k) + g(t[j], x[k]):
Updator := (k, j) -> solve(eq(k,j), yF[k]):

StationaryMatrix := Matrix(n+1, n+1, (k, kk) -> coeff(Updator(k, j), y[2*k-kk])):
SourceVector     := Vector[column](n+1, k -> g(t[j], x[k-1])):
Vars             := Vector[column](n+1, k -> y[k-1]):

StationaryMatrix[1, 1]     := 1:
StationaryMatrix[1, 2]     := 0:
StationaryMatrix[n+1, n  ] := 0:
StationaryMatrix[n+1, n+1] := 1:
SourceVector[  1] := 0:
SourceVector[n+1] := 0:

x := h *~ [$0..n]:
t := tau *~ [$0..m]:

Y0 := [seq(f(x[k]), k=1..n+1)]:
Y := Matrix(n+1, 1, [mu[1](t[1]), evalf(Y0[2..-2])[], mu[2](t[1])]):

for j from 2 to m do
Y[.., j-1];
  Y         := < Y | StationaryMatrix . Y[.., j-1] + SourceVector >:
  Y[  1, j] := mu[1](t[j]):
  Y[n+1, j] := mu[2](t[j]):
end do:

 

plots:-matrixplot(< Vector[column](n+1, Y0) | Y >, style=surface, labels=['x', 't', 'Y'])

 

 


 

Download HeatEquation_Corrected.mw

@tomleslie 

Thank you very much

PS: "... but does in worksheet (honest!)" ... I take your word for it.

@dharr @acer

I didn't pay attention to the "Explore example worksheet" help page.
Thanks for remind me how blind I can be

PS1: I tried something  similar but with using  same name a for the parameters of dsolve AND for the parameters of explore... which returned an error.
Why is the use of different names for the parameters of dsolve and explore necessary?

PS2: I vote up

@Simon45 

By the way, think to check if the condition CFL < 1/2 (strictly less than) is verified. If not increase the spatial step dx or decrease the time step dt.
CFL = lambda*dt/(dx^2) where lambda is the diffusivity cooefficient. Think also that lambda = D/(rho*Cp) for a general heat equation (D = diffusion coefficient)
In your case lambda=1.

Second point: as your scheme is explicit in time (meaning the diffusion term is evaluated at the previous time), you don't have any linear system to solve. You just have to do "time marching".
More precisely, if Y[theta] is the vector which representes the spatial solution atsome time theta, the solution Y[t*dt] given the solution Y[t] is just Y[t+dt]=A*Y[t]+B.
Here A is "your famous tridiagonal matrix", B is the source term g(t, x).
Let M the matrix wich represents the discretized diff(y(x,t), x, x) operator ; then A = dt*I+D where I is the identity matrix.
M is  tridiagonal as soon as one uses (as you did) a centered 3-point finite difference scheme to aproximate diff(y(x,t), x, x) .
M is symmetric if dx is constant over the whole spatial domain.

Third point, an implicit time scheme would result in this expression to update Y:  A*Y[t+dt]=dt*I*Y[t]+B.
In this case Y[t+dt] is obtained by solving a linear system of equations (A is the same A as in the time-explicit case).

An implicit scheme suffers no CFL condition and is unconditionnally stable (meaning every small perurbation of the initial condition is damped as time increases). This means that even a single timestep is enough to obtain the solution at tour final (t=3) physical time.
But be very careful: stability is not enough for the numerical scheme to capture the transient behaviour. 

Last point, updating Y in an explicit schemes is very cheap (matrix vector product) but is submitted to a CFL condition which forces dt to be small.
Doing the same thing for an implicit scheme is costly (solving a linear system) but not restricted to any CFL condition, then you can use larger time steps than with explicit schemes.
There exist better explicit schemes than the simple one you used (less strong CFLcondition).
If you are interested by this topic you will find some informations here:

https://hal.archives-ouvertes.fr/hal-01401125/document  (in English)
 

@Rouben Rostamian  

And NO, I'm not "making this more complicated than it deserves".

I explicitely wrote in my answer to simon45 that I tried to be pedagogical, then your critic is of no matter.

In addition, do I need to remind you that before I answer simon45 you even claimed that "Once you have taken care of satisfying the CFL condition, you will find out that the solution does not deal with tridiagonal matrices at all. ".
So maybe you could find some interest in reading carefully my unnecesary complicated stuff.

@Rouben Rostamian  

Wrong conclusion, you pointed yourself that the CFL conditions has to be verified to obtain a stable scheme:


 

restart;

with(plots):

interface(rtablesize=10):

f:=unapply(-x^2+1,x):

mu[1]:= t -> 0:

mu[2]:= t -> 0:

g    :=(t, x) -> 0:

l   :=  2: T := 3:

n   := 50: m := n*100:

h   := l/n:

tau := T/m:

CFL := tau/h^2;
if CFL >= 1/2 then WARNING("The CFL condition is not verified, the scheme is unstable") end if;

3/8

(1)

FD_T  := k -> (yF[k]-y[k]) / tau:
FD_X  := k -> (y[k+1]-y[k]) / h:
CDD_X := k -> simplify((FD_X(k,j)-FD_X(k-1,j)) / h):

eq := (k,j) -> FD_T(k) = CDD_X(k) + g(t[j], x[k]):
Updator := (k, j) -> solve(eq(k,j), yF[k]):

StationaryMatrix := Matrix(n+1, n+1, (k, kk) -> coeff(Updator(k, j), y[2*k-kk])):
SourceVector     := Vector[column](n+1, k -> g(t[j], x[k-1])):
Vars             := Vector[column](n+1, k -> y[k-1]):

StationaryMatrix[1, 1]      := 1:
StationaryMatrix[1, 2..n+1] := 0:
SourceVector[1]             := mu[1](t[j]):

StationaryMatrix[n+1, 1..n] := 0:
StationaryMatrix[n+1, n+1]  := 1:
SourceVector[n+1]           := mu[2](t[j]):

x := h *~ [$0..n]:
t := tau *~ [$0..m]:

Y0 := [seq(f(x[k]), k=1..n+1)]:
Y := Matrix(n+1, evalf(Y0)):

for j from 2 to m do
Y[.., j-1];
  Y := < Y | StationaryMatrix . Y[.., j-1] + SourceVector >
end do:


plots:-matrixplot(Y, style=surface, labels=['x', 't', 'Y'])

 

 


 

Download HeatEquation_2.mw

First 126 127 128 129 130 131 132 Last Page 128 of 154