Beefy Boxes and Bandwidth Generously Provided by pair Networks
The stupid question is the question not asked

PDEs oo

by Lucero (Novice)
on Nov 02, 2021 at 15:11 UTC ( #11138333=perlquestion: print w/replies, xml ) Need Help??

Lucero has asked for the wisdom of the Perl Monks concerning the following question:

Hello. I'm taking a class, Computational Physics, and I have to solve a few problems with PDEs. The topic was introduced and explained using a book oriented to Python; one of the examples solves the heat equation for a 1-D space and time by leapfrogging (time-stepping) the initial conditions forward in time. I don't have idea how to "translate" that code to Perl and I have searched for similar problems and their implementation in Perl but I found nothing. I'd like to know if exists a module for this kind of problems or website or book where I can find information about this.

Replies are listed 'Best First'.
Re: PDEs oo
by Bod (Curate) on Nov 03, 2021 at 00:26 UTC

    Math::RungeKutta references a book on Partial Differential Equations as one of its sources.
    I've never used it so can't vouch for it but it might provide a starting point.

Re: PDEs oo
by LanX (Sage) on Nov 03, 2021 at 11:04 UTC
    There are no general methods to solve PDEs exactly.

    Most of the time you apply heuristics and/or chose an efficient numerical algorithm (= try-and-error approximation in plain language)

    Numerical approximations mean a lot of number crunching, which are hardly done in Python or Perl.

    The book is most likely addressing external libraries in C or Fortran which interface with Python.

    The heuristic (aka symbolic calculation) parts OTOH can be translated to Perl, if implemented in pure Python.

    IIRC this kind of stuff used to be done in MATLAB or similar tools.

    Cheers Rolf
    (addicted to the Perl Programming Language :)
    Wikisyntax for the Monastery

      > MATLAB -

      I feel more strongly about this than most "GNU" replacements because I really, really don't like the software environment SCATLAB uses; so be sure to mention Octave in the same breath - or maybe as the only - option when it comes to CRAPLAB kind of things :)

        > so be sure to mention Octave in the same breath

        I'm really sorry, but this is still far beyond my wokeness level.


        Cheers Rolf
        (addicted to the Perl Programming Language :)
        Wikisyntax for the Monastery

Re: PDEs oo
by perlfan (Vicar) on Nov 03, 2021 at 06:34 UTC
    Do you want to share the code?

    Here's the 2d version, in C/OpenMP and wrapped in Perl;

    If you want to do this in pure Perl, this can be done rather easily. You're computing a "stencil" (pattern) over space, from the previous timestep. You iterate timesteps until the points in space cease to change by some difference (epsilon). You can also do this in PDL. Other relevant terms: Jacobi, Gauss-Seidel, SOR.

      from numpy import * import matplotlib.pylab as p from mpl_toolkits.mplot3d import Axes3D Nx = 101; Nt = 3000; Dx = 0.03 ; Dt = 0.9 KAPPA = 210. ; SPH = 900. ; RHO = 2700. #conductivity,specfheat ,densi +ty T = zeros( (Nx , 2) , float ); Tpl =zeros( (Nx , 31) , float ) print ("Working, wait for figure after count to 10") for ix in range ( 1 , Nx - 1) : T[ ix , 0] = 100.0 ; # initial temper +ature T [ 0 , 0 ] = 0.0 ; T [ 0 , 1 ] = 0. #first and last points a t0 T[Nx - 1 ,0] = 0. ; T[Nx - 1 ,1] = 0.0 cons = KAPPA / ( SPH*RHO) *Dt / ( Dx*Dx) ; # constant m=1 # counter for rows , one every 300 time steps for t in range ( 1 , Nt ) : # time iteration for ix in range ( 1 , Nx - 1) : #Finitedifferences T[ ix , 1] = T[ ix , 0] + cons *(T[ ix + 1 , 0] + T[ ix - 1, 0] + - 2.*T[ ix , 0 ] ) if t %300 == 0 or t == 1: #for t = 1 and every 300 time steps for ix in range ( 1 , Nx - 1 , 2) : Tpl [ ix , m] = T[ ix , 1] print (m) m=m+1 # increasemevery 300 time steps for ix in range ( 1 , Nx - 1): T[ ix , 0] = T[ ix , 1]# 100 p o s i + t o n s a t t =m x = list ( range ( 1 , Nx - 1 , 2) ) # p l o t every o t h e r x p o i + n t y = list ( range ( 1 , 30) ) # every 10 p o i n t s i n y ( time ) X, Y = p.meshgrid ( x , y ) # g r i d f o r p o s i t i o n and time def functz ( Tpl ) : # Function r e t u r n s t e m p e r a t u r e z = Tpl [X, Y] return z Z = functz ( Tpl ) fig = p.figure( ) # c r e a t e f i g u r e ax = Axes3D ( fig ) # p l o t s a x i s ax.plot_wireframe(X, Y, Z , color = 'r') # red wireframe ax.set_xlabel('Position') # l a b e l axes ax.set_ylabel('time') ax.set_zlabel('Temperature') ) # shows f i g u r e , c l o s e Python s h e l l print ("finished")
      Switching to Perl, I can change the matrix elements but I dont know how to change the dimensions of the matrix.
Re: PDEs oo
by chafelix (Novice) on Nov 02, 2021 at 19:39 UTC
    There is a very good reason why such work is done in FORTRAN: Neither Python, nor Perl are the right tools for PDEs. The reason is that PDEs can be expensive to solve and you want to choose the RIGHT code (e.g. implicit schemes for stiff problems) that runs FAST. Where Perl or any other interpreted language for that matter has a place is as a front-end or backend, for instance a graphical interface to set the parameters and possibly parse the results and create say an Xmgrace graphics file. So the answer is I cannot be certain noone has done a perl implementation of a pde solver, but honestly I'd consider this rather a waste of time.

      While not necessarily a "wrong" answer, I think this is kind of a bad answer to the original question.

      If you're doing this "seriously" for something, then yes you're going to want to throw something bare metal and these days probably GPU powered at it. However for the purposes of learning how they work and how they're implemented and solved I think something like perl (or python; bleh) is going to be perfectly reasonable and performant for the type of examples you'd encounter in problems. People used to solve these things with slide rules and pencil/paper. If your (again, experimental learning) implementation to solve the hairiest problem in your textbook takes 2 minutes rather than 10 seconds that's not really going to be an issue the tens of times you might be running it.

      And as to the OP not being able to understand the python examples: if you're going to be programming seriously go ahead and take the time to get a reading familiarity of python (the Pascal of the naught-ies . . .). Even if you stick with Perl as your primary language it's worthwhile to be able to digest things in other languages because there's going to be times like this where you can't find something off the shelf. To say nothing of being exposed to other language approaches helps you improve your programming in general.

      The cake is a lie.
      The cake is a lie.
      The cake is a lie.

        That's what Im doing; If I don't know how, or can't find a way, to do it in Perl, I switch to another language.Not the first time, it happens. At the moment, I think I found an answer to my problem but I'm still working on it.

Log In?

What's my password?
Create A New User
Domain Nodelet?
Node Status?
node history
Node Type: perlquestion [id://11138333]
Approved by erzuuli
and the web crawler heard nothing...

How do I use this? | Other CB clients
Other Users?
Others lurking in the Monastery: (3)
As of 2022-08-12 23:26 GMT
Find Nodes?
    Voting Booth?

    No recent polls found