mmcdara

6149 Reputation

17 Badges

9 years, 73 days

MaplePrimes Activity


These are replies submitted by mmcdara

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

@Rouben Rostamian  

You did not carefully read my explanations and you make a confusion between the mathematical formalism of the problem and its practical solution.

Mathematically a PDE problem with bcs in a bounded domain is always made of two things : a PDE that describe what happens in the OPEN domain, plus bc that describe what happens on its BOUNDARY.
From this point of view you're right saying the matrix is tridiagonal.

Now look to any book concerning "practical implementations" of finite differences, finite volumes or finite elements methods.
For instance (apologies, I do not know an equivalent reference in english but I'm sure there are many)
https://ebook.chapitre.com/ebooks/methode-des-elements-finis-9782746296664_9782746296664_2.html  (in French)

One of the question is: how to treat the boundary conditions when a global matrix system is  build by assembling together, in a very systematic way, small local matrices ?
One of the simplest way to proceed is to build a matrix in the CLOSED domain, as if no boundaries existed, and after that to correct this matrix by "adding" the boundary conditions. The trick is to introduce fictitious points outside of the domain in order to extend the discretization of the spatial differential operators up to the boundary, and next to use the boundary conditions to eliminate these fictitious points.

To fix the ideas, let's suppose that the spatial domain is discretized at point x[0], ...x[N+1].
The "interior" matrix has size N by N and you have 2 extra relations describing bc at points x[0] and x[N+1].
From the mathematical point of view there is no need to try to solve the problem on the boundaries (particularly if you have set Dirichlet conditions), but from practical point of you it's better to deal with a (N+2) by (N+2) matrix which includes the boundary conditions.

This difference between the mathematical and the practical point of view is exactly the same we met in spline approximation: would you say the matrix of a cubic spline is tridiagonal or not?
If you do not consider the first and last points, yes it is, but no longer as soon as you include them in the global matrix.
 

First 97 98 99 100 101 102 103 Last Page 99 of 125