mmcdara

7891 Reputation

22 Badges

9 years, 56 days

MaplePrimes Activity


These are replies submitted by mmcdara

@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.
 

@Rouben Rostamian  

Hi,

I agre that the stability of the FTCS (Forward Tme, Centered Space) scheme for parabolic equations requires a CFL conditions is satisfied.
As you pointed, this condition (lambda*dt/dx^2 < 1/2; dt = time tep, dx= space step, lambda [=1 for the OP] is the diffusivity coefficient) is not verified in this problem (dt/dx^2 = 15/2).

So the choice of dt, for instance, is very poor

On the other hand, I do not agree with you writting "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."

I think that there is a general misunderstanding around this "tridiagonal" term (maybe the OP was a little bit confused or at least did not expressed that clearly). The tridiagonal matrix the OP refers to is the matrix of the "updator" used at each time step.
The structure of this matrix  is quasi tridiagonal. It corresponds to a tridiagonal N by N matrix,  with 1 colum of zeroes added to the left and to the right, with a row of zeroes added to the bottom and the top, and the elements [1, 1] and [N+2, N+2] set to 1 (to account for the boundary conditions).
But, in the OPEN spatial domain, the system is tridiagonal because a Centered finite difference approximation of diff(y(x,t), x, x) is used: this is true whatever the value of the CFL.

You wrote "The tridiagonal system enters the picture when you do the time discretization through the implicit (or backward) finite differences": not exactly

  • with an explicit scheme a stationnary one writes that the value of y(x, t+dt) =y(x, t)+ something (the moral drawback of this simplicity is the respect of the CFL condition)
    More precisely, et Y[t] the vector column < seq(y(x[k], t+dt), k=0..n) >, then Y[t+dt]=A*Y[t]+B
    A represents the centered finite difference approximation of the diffusion operator and B extra source terms.
    A is a quasi tridiagonal matrix
     
  • with an explicit scheme, one writes A*Y[t+dt]=C, where C is source term vector (not the B as before).
    A is a quasi tridiagonal matrix withs exactly the same shape as the A used in the explicit scheme.
    Note that Y[t+dt] is nox obtained by solving a linear systeme as it was previously obtained by doing a simple matrix-vector product

So the main difference lies essentially in the way A is used, but whatever the type of scheme (explicit or implicit), A is always a quasi tridiagonal matris

 

@Rouben Rostamian  

Thank you for giving a solution and the reference that inspired it when it is not entirely your own.
Doesn't seem so common here.

@ebrahimina 

I'm definitively not the most competent to answer this question.

I  suppose you're looking for online sources, not books to buy?
If so, you could try this one 

https://www.maplesoft.com/documentation_center/maple18/ProgrammingGuide.pdf

Because  I have the paper manual delivered by my old Maple 8 and written by the same authors (more or less), I suppose this update still remains a very good reference.

 

@Mariusz Iwaniuk 

For information, it worked with Mape 2015=

 

restart:

interface(version)

`Standard Worksheet Interface, Maple 2015.2, Mac OS X, December 21 2015 Build ID 1097895`

(1)

solve(0.1 = 23.714*(-0.93205)^2/(20.3+61.4*.884^x), x)

-8.976314317

(2)

 


 

Download solve2015.mw

@vv 

The OP wrote " Let the t-axis go from 0 to 2*Pi.":this also could be interpreted as "I want to restrict the plot to 0..2*Pi".

But I don't want to give the impression that I'm right at all costs or to be polemic, so I'm giving up this thread.

See you on another one :-)

@vv 

SORRY, 

I used add, verified it worked (returned infinity), saved (or not?), changed add into sum, did probably some orther thing 

and finally screwed up everything.

Of course it's add

@vv 

 

f:=t -> t^2;

a := n -> 1/Pi*int(cos(n*t)*f(t),t=-Pi..Pi):

b := n -> 1/Pi*int(sin(n*t)*f(t),t=-Pi..Pi):

S := (n,t) -> a(0)/2+add(a(k)*cos(k*t)+b(k)*sin(k*t),k=1..n):

plot([f(t),S(8,t)], t=-Pi..Pi);

 

proc (t) options operator, arrow; t^2 end proc

 

 

 


 

Download WhyNot.mw

@Carl Love   @vv

The initial question was about the Fourier Coefficients (FC for short) of the function f defined by f:t -> t^2.
This function not being L1 nor periodic these FC do no not exist.

The first question then becomes: "How to define a periodic L1 function g such that g(t) = f(t) over some real bounded interval [a, b]  (then g is L1 over [a, b]) and that g replicates periodically outside [a, b] ?
You can choose [a, b] =[0, 2*Pi] or anything else, for instance [a, b]=[-Pi, +Pi]... but your choice impacts the FC (in the latter case the function g being symetric around 0, all sine FC are null).
Before calculating the FC one must define what this function g is.

The second point is that restricting f to [a, b] is equivalent to multiply f(t) by G(t) where G is the "gate function" defined by
Heaviside(t-a)-Heaviside(t-b) whose Fourier transform is, after normalization, a cardinal sine (~sin(t)/t)
The Fourier Transform of f.G is TF(f)*TF(G)  where * denotes the convolution product and TF the Fourier Transform.
Even if the spectrum of f was bounded, the spectrum of f.G. would not be bounded because of this cardinal sine. 
Because of the truncation of the signal its spectrum becomes bounded. Its sampling cannot verify rigously Shannon's sampling theorem and aliasing appears.
Because g is a periodic, so is its spectrum, and the "problem" introduced by the convolution with TF(G) periodically replicates, making "noisy" the extreme part of the spectrum.

All the art in signal processing is to define a good observation range [0, T], a good sampling frequency, a good windowing to limit aliasing, and a good [a, b] domain.
One choice is [a, b]=[0, T], another is [a, b]=[0, 2T] with g(t)=f(t) if t is in [0, T] and g(t)=0 if t is in [T, 2T].

The Gibbs phenomenon you talked about is more related to the discontinuities of g than to the points discussed above.

Here is the first reference I found on the web concerning discrete fourier transform and aliasing. I do not say it's the best one, but more detailed analysis can be found in any book concerning Discrete Fourier Transform (which is finally the practical application to Fourier's series expansion to true signals).

https://www.tcd.ie/Physics/research/groups/mobius/teaching/Numerical%20methods%20III/PY4C01_2012_L5.pdf

This is more specifically for vv.
I agree that "The standard periodic extension of a function f defined in [0,T) is simply x |--> f(x - floor(x/T)*T)." but this formula is impractical to understand what happens when you compute Fourier coefficients or do FT, FFT or DFT.
The correct framework to do the analysis and to understand why and how periodization affects the FC is to use "Dirac comb" (not a function but, formally, a distribution)

See for instance (the term "sha" is replaced by "Dirac comb")

https://en.wikipedia.org/wiki/Dirac_comb
https://pdfs.semanticscholar.org/5d3f/61dfe1cfbf59bd25d6bff343f0cc520bd955.pdf

 

 

@Rouben Rostamian  

Agree: it's not clear if the OP is "asking for help with the mathematics of the problem, or how to solve it in Maple"/

To answer your question "So what's the problem with viewing the function f(t)=t^2 as periodic?": perization of f(t), assuming t is in [0, T] is not trivial.
At least two classical methods are used based on the "sha function" ( an infinite sequence of translated Dirc distributions delta(t-kT) where k is in Z:

  1. form the convolution product of f(t) by sha[T](t)  (here sha[T](t) is the "sha function" of period T)
  2. define g(t) = f(t) if t is in [0, T] and 0 if t is in [-T, 0] and form the product of g(t) by sha[2T](t)

In the first way there can be a strong "spectrum folding" (aliasing effect)  which biases the estimation of thr Fourier coefficients. The slower the decreasing of the power spectrum, the higher the bias.
Windowing (in t space ot frequency space) is a classical method to limit aliasing.

The second approach also reduces aliasing in a slightly simpler way.

At the end the true problem is: what is the more reasonable way to periorize f given the problem at hand?

Hi,

It works if you describe the summation differently
 

restart:

interface(version)

`Standard Worksheet Interface, Maple 2015.2, Mac OS X, December 21 2015 Build ID 1097895`

(1)

L:=[1,2,0,4]:

try    # bug

  s := sum(1/u, u in L);

  catch:  s:=infinity

end try;

 

infinity

(2)

 


 

Download sum.mw

@Carl Love 

Tank you Carl.

If it's not too much to ask, why geometry:-triangle(...) returns an error? it seems that geometry Toto does not behave like the "ordinary" packages

@acer 

I  just tested your code with my wife's Maple 2016 version: I agree vv, it's very nice.

I will be sure to let you know if any adjustments (which I could not have made on my own) could indeed be given.

Thank you for all your involvement

PS: in response to a previous mail of yours where you wrote "On the help page for DocumentTools:-Tabulate the first sentence in the paragraph where it describes the coloroption says, "The foreground color for text." ..., you're probably right writting " I interpret "text" as meaning entries which are strings and not expressions. Which means that the behavior is just as documented. Not a bug " , I didn't payed enough attention to this point.

@tomleslie 

You didn't missed anything.

Thanks for the code.

@acer 

Thank you acer.
I have only Maple 2015 right now and I can't test, but the display suits me perfectly.

PS: In case you want to understand my motivations about   his question, my answer to Tomleslie will give you a little idea.

Thanks again

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