patchGNU Octave - Patches: patch #7929, Partial Differential Equation...


patch #7929: Partial Differential Equation Solver in two dimensions

Submitted by:  None
Submitted on:  Sat 19 Jan 2013 12:19:30 AM UTC  
Category:  None Priority:  5 - Normal
Status:  In Progress Privacy:  Public
Assigned to:  Carlo de Falco <cdf> Originator Email:  -email is unavailable-
Open/Closed:  Closed

Add a New Comment (Rich Markup)

You are not logged in

Please log in, so followups can be emailed to you.


( Jump to the original submission)

Fri 29 Mar 2013 10:22:35 AM UTC, comment #18: 


As Jordi seems to be too buisy to review this patch I'll review it during the Easter weekend and remove the item n the GSoC 2013 list of ideas.

Would you also be interested in implementing the BVP*C functions?


Carlo de Falco <cdf>
Project MemberIn charge of this item.
Thu 14 Mar 2013 09:56:22 PM UTC, comment #17: 

I see that you published the SoC ideas for 2013.  My code completes the workitem:

"Implement solver for initial-boundary value problems for parabolic-elliptic PDEs in 1D"

Performance is a bit of concern, but that should be taken care of by the JIT compiler, I think.

Mark Wistrom <mwistrom>
Wed 13 Mar 2013 06:31:24 PM UTC, comment #16: 

Yes, the code in the .tar file satisfies the requirements. The functionality is very similar to pdepe in Matlab. This version should run the same code as Matlab's with little or no modification and produce identical or better results.

Mark Wistrom <mwistrom>
Wed 13 Mar 2013 03:41:07 PM UTC, comment #15: 
Jordi Gutiérrez Hermoso <jordigh>
Project Administrator
Thu 14 Feb 2013 07:59:25 PM UTC, comment #14: 

I'll add this comment to the patch so that we do not lose the information.

Let me know if you have any more comments on the code.


-----Original Message-----
From: Dongxu Zhou
Sent: Thursday, February 14, 2013 11:34 AM
To: Mark Wistrom
Subject: RE: [patch #7929] Partial Differential Equation Solver in two dimensions

Hi Mark,

It turns out that I mistakenly transposed u without noticing it.  So the pdepe function is correct. 

Thank you for developing the code!

Mark Wistrom <mwistrom>
Thu 14 Feb 2013 06:23:11 PM UTC, comment #13: 

Let me check it out tonight.

THANK YOU for the comment.  My goal is to make this a very useful feature.


Mark Wistrom <mwistrom>
Thu 14 Feb 2013 05:56:16 PM UTC, comment #12: 

Hi Mark,

Thank you for the code. This is an function I am desperately looking for.

While running the code, I noticed that there might be a bug: if I provide different n values for x and t,

  x = linspace (0, 1, 10);
  t = linspace (0, 1, 20);

the surf function won't run because of dimension mismatch.  Is there a matrix transpose problem for the code?


Dongxu <zhoudx>
Tue 12 Feb 2013 06:35:40 AM UTC, comment #11: 

attached files:

pdepe.m     solve PDE
pdeval.m    evaluate solution at different points

osb##.m     examples

(file #27448)

Mark Wistrom <mwistrom>
Tue 12 Feb 2013 03:35:01 AM UTC, comment #10: 

The '1' at the top of the files is due to incompetence. The '1's have all been removed.

My call is core function, but I won't argue hard for that. Matlab has this as a core function. It seems to me that pdepe is a "nice to have" or "handy" function.  I get the impression that it will be used for education, prototyping or "single shot" calculations. Not something that would be used in a heavily used or "industrial" scenario, if that makes any sense. Ease of discovery and ease of use seem to be important. But I don't know how you are positioning Octave.
I think it will probably used mostly for education. See example 5: "Cooking the Turkey."

I've read the style issues you pointed me to and made the following changes:

1. change 'inline' to @
2. change all "ends" to "endif", etc
3. change comments to all "%"
4. add space between function and argument to distinguish functions from indexing
5. use single single quoted strings
6. pare lines to under 80 char (except for the function definition in the 'help'.  I think I got it right.)
7. use !=
8. use spaces, no tabs
9. use Unix newlines
10. put spaces after ','
11. use lower case names for functions and variables
12. enclose "if" and "switch" test conditions with ()

I have removed much of the duplicate code. There was more than I remembered.  :)

So the trivial loops are in the examples, which I do know what to do with. I have not changed the examples to use more efficient coding, because I don't know where those examples should go. I can change them, if needed. Do they get used as tests? Do I provide tests with the code? Also, the examples return a single value for the 'c,f,s' parameters, not vectors, just like Matlab's functionality. This seems to have a pretty big performance hit, which I discussed below.

I'll post the functions and examples soon.


Mark Wistrom <mwistrom>
Thu 07 Feb 2013 12:42:06 AM UTC, comment #9: 

Thank you for your time looking this over. It is excellent feedback. I will address most if not all of your comments and concerns. First a couple of overall comments.

1. The mathematics of the implementation. When I first looked at the problem, I thought it would be easy if not trivial. The more I looked into it, the more I realized why this had not been done before. The paper it is based on takes some time and effort to get up to speed and is a little opaque. I think I have the proper mathematics in place based on the testing that I have done. It agrees with Matlab in most cases, and in the Singularity case the error "looks" right. i.e. it scales nicely as n^-2, is smooth with no large jumps, and is better than Matlab's.

2. The "trivial remarks" that you have made are valid and probably not trivial. I will work on them.

3. Code duplication. Yes, there is code duplication. I put that in there as part of prototyping the functionality to ensure I was making head way. I'll address that.

4. Performance. It is slow. I know. I used the profiler (That thing is pretty nice BTW), and the loop where the user function is called to build the vectors is causing the problem, which is what you would expect. I see two issues: 1. The call into the user function is per each finite element at each time step. If the user function were passed in the vector of values for x,u,dudx and returned a vector, that speeds things up at the expense of being different from Matlab's implementation.  2. The loop itself could be compiled.

5. Matlab's pdepe is said to solved parabolic AND elliptic equations (at least one needs to be parabolic). I haven't tested the elliptic part of the my code. My code sets up the equations and then kicks off daspk (Matlab uses ode15s) to do the time integration. So I'm not really sure what to do here. Document things most likely.

I'll make code changes and post back with specific answers.

Thanks again,

Mark Wistrom <mwistrom>
Wed 06 Feb 2013 09:50:56 PM UTC, comment #8: 

Sorry, I messed up the code fix, I meant

function [u] = osb22ASol(x , t)
  P = 1.1;
  C = ifelse (x < 0, 0.1, 1);
  u = log (C.*x(:) + t(:)' + P)

Jordi Gutiérrez Hermoso <jordigh>
Project Administrator
Wed 06 Feb 2013 09:49:35 PM UTC, comment #7: 

Okay, so I'm starting to think about importing your code into Octave...

I won't comment on the mathematical aspects of it, so please be patient with the apparently trivial remarks I have to make about its superficial appearance.

Why are you putting 1; at the top of files, in effect using script files and no function files?

Should this be a package or a core function?

There are many stylistic issues to also discuss. We prefer a different Octave house style to Matlab style. Read some of Octave's source code to glean this style and this:

I see the following problems with your code:

  • Using inline("foo(a,b)", "a", "a") instead of @(a,b) foo(a,b)
  • No explicit end statements (endfor, endif, endwhile...)
  • % for comments instead of ##
  • No space between function names and opening bracket: sin (x), not sin(x) (but note, no space for indexing, just a visual cue to distinguish function calls from indexing).
  • Single-quoted strings. Except in some cases, we generally prefer double-quoted strings.
  • Exceedingly long lines. Lines should not be more than 80 characters wide. Note that unlike Matlab, Octave does not require continuation characters if there's an opening bracket that hasn't been closed on this line yet. Use this feature to your advantage when breaking up lines.
  • Using ~= instead of the more common !=
  • No spaces after commas.
  • Inconsistent use of tabs and indentation. Just don't use any tabs at all and use two-space indentation.
  • Windows newlines (/r/n). We use Unix newlines (just /n).

Here is a more profound problem: I see a lot of code duplication here. Code duplication is very bad. In the future as we need to maintain this code, code duplication means that we have to hunt down all instances of the same code and patch it in the same way. Please move as much common code as you can into common functions shared by others. If necessary, use private directories and or subfunctions.

Another problem, I see some very trivial loops that can easily be replaced. Instead of

function [u] = osb22ASol( x , t)
    C1 = 0.1;
    C2 = 1.0;
    P = 1.1;
    for i=1:length(t)
        for j=1:length(x)
            if x(j)<0
                u(i,j) = log( C1 * x(j) + t(i) + P);
                u(i,j) = log( C2 * x(j) + t(i) + P);

this could easily be

function [u] = osb22ASol(x , t)
  P = 1.1;
  C = ifelse (x < 0, 0.1, 0.2);
  u = log (C*x(:) + t(:)' + P)

This is just from a quick glance of the code, sorry I didn't delve deeper into the actual mathematics of it.

Jordi Gutiérrez Hermoso <jordigh>
Project Administrator
Tue 29 Jan 2013 10:04:23 PM UTC, comment #6: 

Examples SB21 and SB27 are both working properly now that I have found and corrected the errors in the paper.

Same thing with all the equations where there is a singularity though. Matlab's error is about an order of magnitude larger, isn't smooth and doesn't drop nicely as the number of finite elements increases.

I think they didn't implement what is in the paper.  <shrugs>

Mark Wistrom <mwistrom>
Tue 22 Jan 2013 08:37:58 PM UTC, comment #5: 

User error.

Those were supposed to be example files for Matlab to show the results of the same examples.

I've attached the file again with the correct examples with a different, less ominous filename.


(file #27313)

Mark Wistrom <mwistrom>
Tue 22 Jan 2013 08:10:21 PM UTC, comment #4: 

matlabfiles.tar is empty. What was supposed to be in it? I hope you are not trying to distribute Matlab source code without permission.

Jordi Gutiérrez Hermoso <jordigh>
Project Administrator
Tue 22 Jan 2013 01:51:55 AM UTC, comment #3: 

attached the latest pdepeO.m file and the Matlab and SB examples.

(you all probably know this :), but... )
usage in octave for right now:


(file #27308, file #27309)

Mark Wistrom <mwistrom>
Tue 22 Jan 2013 01:34:49 AM UTC, comment #2: 

The five pdex# examples for Matlab's pdepe function work for my code and seem to produce the same results. The 10 examples from the SB paper produced mixed results when compared with the exact analytical solution which were provided. The biggest surprise is that the error from MATLAB's implementation (tested on R2011a) is about an order of magnitude worse than my code for the examples where there is a singularity. This is interesting because the entire purpose of the paper was to provide a general PDE solver for the cases where there is a singularity. I'm not sure what MATLAB did, but their errors are clearly larger and they have some irregular jumps/features in the error near the origin that are absent from my code. Perhaps this from user error or something else I did, but I am at a loss to explain or fix it.

The examples are numbered by their subsection from the paper. The comparisons between the two implementations are below:

SB11  m=0 
Both produce the same result

SB21  m=2
SB25  m=1
Matlab says that the index is too big. daspk returns without being able to find initial conditions.  No solution.

SB22  m=2
SB23  m=2
SB26  m=1
SB28  m=1
In these four examples Matlab produces a solution that is about an order of magnitude worse than my solution and has strange features around the origin.  NOTE: There is a typo in the analytical solution for SB23, 4xt should be 6xt.

SB24  m=2
Solvers generally give the exact solution for this problem and my code does. Matlab states that the index is too big.

SB27  m=2
my code and Matlab produce the same strange (incorrect) solution. I do not know why.

SB29  m=1
Mathematica, Matlab and my code produce the same result which is quite different from the exact solution given in SB. I checked that the exact solution solves the equation.  I do not know what is going on here.

Mark Wistrom <mwistrom>
Sat 19 Jan 2013 12:24:01 AM UTC, comment #1: 

Although it is in two dimensions, I wish the title were better now that I have read it.  Also, I'm logged in now.  This was my original submission.

Mark Wistrom <mwistrom>
Sat 19 Jan 2013 12:19:30 AM UTC, original submission:  

This is an equivalent implementation of MATLAB's pdepe PDE solver based upon a paper by Skeel and Berzins. The solver has two basic parts.

1. configuring the system based on the coordinates of the system and the presence of a singularity at the origin of the polar and spherical problems.

2. providing the function for integration by an ODE/DAE solver.

I have several more examples from both the examples from Matlab and from the SB paper.

I'm at the point where I think I'm winding down on implementation and looking for a list of things to do to it to completely finish.



(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:

No files currently attached


Depends on the following items: None found

Items that depend on this one: None found


Carbon-Copy List
  • -email is unavailable- added by cdf (Posted a comment)
  • -email is unavailable- added by zhoudx (Posted a comment)
  • -email is unavailable- added by jordigh (Posted a comment)
  • -email is unavailable- added by mwistrom (Posted a comment)
  • -email is unavailable- added by None (Submitted the item)

    Do you think this task is very important?
    If so, you can add your encouragement to it.
    This task has 0 encouragements so far.

    Only logged-in users can vote.




    Follow 19 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2013-03-31 cdf Open/ClosedOpen => Closed
    2013-03-31 cdf Attached File#27302 => Removed
    2013-03-31 cdf Attached File#27301 => Removed
    2013-03-31 cdf Attached File#27300 => Removed
    2013-03-31 cdf Attached File#27299 => Removed
    2013-03-31 cdf Attached File#27309 => Removed
    2013-03-31 cdf Attached File#27308 => Removed
    2013-03-31 cdf Attached File#27313 => Removed
    2013-03-31 cdf Attached File#27448 => Removed
    2013-03-29 cdf Assigned toNone => cdf
    2013-02-12 mwistrom Attached File- => Added pdepe.tar, #27448
    2013-02-06 jordigh StatusNone => In Progress
    2013-01-22 mwistrom Attached File- => Added ExamplesForMatlab.tar, #27313
    2013-01-22 mwistrom Attached File- => Added Matlabfiles.tar, #27308
        Attached File- => Added ofiles.tar, #27309
    2013-01-19 None Attached File- => Added pdepe.m, #27299
        Attached File- => Added Opdex4.m, #27300
        Attached File- => Added Opdex1.m, #27301
        Attached File- => Added Osb24.m, #27302

    Back to the top

    Powered by Savane 3.5