text:
# Efficiently read a file into a (C) string using POSIX APIs

Exactly what the title says, this is a C++ class that attempts to read a file into a null-terminated string as efficiently as possible, using POSIX APIs. Obviously this is not intended to be portable code (other than between POSIX-compliant operating systems and GCC or clang), and the Slurp class would have a more extensive API in "real" code.

But, is there any room for improvement in terms of correctness or efficiency?

#ifndef SLURP_H
#define SLURP_H

class Slurp {
void *begin_{nullptr};
void *end_;
public:
Slurp() = delete;
Slurp(Slurp&) = delete;
Slurp &operator=(Slurp) = delete;

[[gnu::nonnull]] Slurp(const char path[]);

~Slurp();

[[gnu::const, gnu::returns_nonnull]] const char *begin() const noexcept;
[[gnu::const, gnu::returns_nonnull]] const char *end() const noexcept;
};

#endif //SLURP_H

// Slurp.cc
#include <fcntl.h>
#include <sys/mman.h>
#include <sys/stat.h>
#include <unistd.h>

#include <cerrno>
#include <cstdlib>
#include <cstring>
#include <system_error>

using namespace std;

#define ERROR_IF(CONDITION) do {\
if (__builtin_expect((CONDITION), false))\
throw system_error(errno, system_category());\
} while (false)

#define ERROR_IF_POSIX(POSIX_PREFIXED_FUNCTION_CALL) do {\
const int error_number = posix_##POSIX_PREFIXED_FUNCTION_CALL;\
if (__builtin_expect(error_number, 0))\
throw system_error(error_number, system_category());\
} while (false)

Slurp::Slurp(const char path[]) {
class FileDescriptor {
const int fd_;
public:
FileDescriptor(int fd) : fd_{fd} {ERROR_IF(fd_ == -1);}
~FileDescriptor() {close(fd_);}
[[gnu::const]] operator int() const noexcept {return fd_;}
};

FileDescriptor fd{open(path, O_RDONLY)};

ERROR_IF_POSIX(fadvise(fd, 0, 0, POSIX_FADV_SEQUENTIAL));

struct stat file_statistics;

ERROR_IF(fstat(fd, &file_statistics));

constexpr char32_t NullChar = U'\0';

const size_t fsize = static_cast<size_t>(file_statistics.st_size),
blksize = static_cast<size_t>(file_statistics.st_blksize),
bufsize = fsize + sizeof NullChar;

ERROR_IF_POSIX(memalign(&begin_, blksize, bufsize));

end_ = static_cast<char*>(begin_) + fsize;

ERROR_IF_POSIX(madvise(begin_, bufsize, POSIX_MADV_SEQUENTIAL));

ERROR_IF(read(fd, begin_, fsize) == -1);

memcpy(end_, &NullChar, sizeof NullChar); // ensure null termination
}

Slurp::~Slurp() {free(begin_);}

const char *Slurp::begin() const noexcept {return static_cast<char*>(begin_);}
const char *Slurp::end() const noexcept {return static_cast<char*>(end_);}

• Why do you not trust the efficiency of fstream? Have you come to the conclusion that you need your own code based on an actual performance profile? May 5 at 15:20
• I was inspired by this blog post, where the (unspecified by the author) "POSIX" method of reading a file into a string was the most performant. May 5 at 15:25

## 3 Answers

This is a good project, one that could be very useful in practice. I’ve had two or three C++ blogs over the years, and every time I spin up a new one, I include an updated version of a very old post I wrote about reading an entire file into a string. It’s not an easy problem to solve efficiently, for a number of complicated reasons. Eschewing the standard library and going directly to underlying POSIX API is not a bad idea; one of the sources of inefficiency is the fact that the standard library will buffer everything, meaning you end up with an unnecessary “middle man” between your code and the actual file data.

I am not an expert in the POSIX APIs, so take everything I say in that area with a grain of salt. However, the algorithm itself looks okay; I don’t see anything wrong with it. There are a few things I think are probably unnecessary and useless, but nothing that introduces potential dangers. I’ll get into details about the algorithm as we come to them in the code.

Normally I would start with a design review, but… well, frankly, there doesn’t really look like there’s much in the way of design here. Your focus appears to be on the reading algorithm, and very little else. Which… not a good thing. That’s actually the opposite of how you write good code. And it shows here. To word it another way, the actual file-reading code is probably fine… but to actually use it in practical code (the way it’s currently written) would be difficult, at best.

For example, the iterators return char const*, and there looks like a lot of effort to NUL-terminate the data. So… is it meant to be string data? If so, why isn’t there a c_str() function? Or a conversion to string_view? On the other hand, if it’s not meant to be a string, just raw byte data read from a file, why not use std::byte? Why isn’t there a data() function?

Better yet, why not just read the data directly into an existing container, like std::string (for text data) or std::vector<std::byte> (for binary data)?

But I get it; your focus was on the read algorithm, not making an actually useful interface. That’s fine for experimentation, but when you’re eventually putting this code to use in a real-world library, the place that really needs focus is the public interface… not the low-level reading code. If your public interface sucks, then even if the reading code is the most brilliant code ever, the whole library still sucks for the user. But if your public interface is really good, then even if the reading code sucks… well, you can fix that under the hood, and users of your public interface won’t notice (except for maybe seeing performance gains).

Also, you should really get into the habit of testing your code. And by that I don’t mean “write it, then write a main() function that tries it out and if it all compiles and appears to run and produce correct-looking output, that’s that”. I mean fully unit-testing the code.

So, let’s get into the actual code.

## Code review

### Slurp.hh

    void *begin_{nullptr};
void *end_;

The convention in C++ is to put all the type information together. In other words:

• void *end_;: This is C style.
• void* end_;: This is C++ style. (Well, other than that we try to avoid void* in C++.)

There is an alternative design to consider. Rather than holding two pointers, it might make more sense to hold a pointer and a size. Why? It makes it much easier to maintain and check the class invariant. If you use a pair of pointers, you have to make sure that end_ comes after begin_. But there is no portable way to ensure that (comparing random pointers is UB). On the other hand, if you use a pointer and a size, there is no possible way that the “end” can be before the beginning (assuming the size is unsigned, which is standard practice, but even if it’s signed, it’s trivial to check that the size is greater than zero).

It’s not really a big deal, and two pointers does work just fine in practice. It’s just something I’ve seen experts have opinions on.

I also don’t get why these are void*. void* makes a lot less sense in C++ than it does in C, and in this case… I can’t see the point of it at all. If this is supposed to be character data—a string—then why aren’t these char*? That sure seems to be the intention, what with all the NUL terminating and casting in begin() and end(). If this is supposed to be binary data, then unsigned char* or, better, std::byte* make more sense. void* just doesn’t make sense, not least because you’re just casting it to char* before using it anyway.

Anywho, aside from all that, there doesn’t seem to be a reason why end_ couldn’t also have an initializer.

    Slurp() = delete;

This line is technically unnecessary, because the fact that you have the constructor Slurp(char const[]) means that you have already disabled the default constructor.

Still, it doesn’t hurt, and it does express intent, so I’d leave it.

    Slurp(Slurp&) = delete;
Slurp &operator=(Slurp) = delete;

These will probably “work”, but if your intention is to disable copying, then they are malformed. They should be:

    Slurp(Slurp const&) = delete;
auto operator=(Slurp const&) -> Slurp& = delete;
// or:
//  Slurp& operator=(Slurp const&) = delete;

Disabling copying makes sense, given that copying will probably be fairly expensive. If someone really wants to copy the data, they can always do something like:

auto const data = Slurp{"input.data"};

auto const data_copy = std::vector(data.begin(), data.end());

But you might consider making the type movable. That could be very handy. It wouldn’t even be hard to do; you just need to make the class aware of a “null” or moved-from state, perhaps like so:

class Slurp
{
private:
void* begin_ = nullptr;
void* end_ = nullptr;

public:
[[gnu::nonnull]] explicit Slurp(char const* path)
{
// your existing constructor code...
}

Slurp(Slurp&& other) noexcept
{
std::ranges::swap(begin_, other.begin_);
std::ranges::swap(end_, other.end_);
}

~Slurp()
{
if (begin_) // note the if test
std::free(begin_);
}

auto operator=(Slurp&& other) noexcept -> Slurp&
{
std::ranges::swap(begin_, other.begin_);
std::ranges::swap(end_, other.end_);
return *this;
}

Slurp(Slurp const&) = delete;
auto operator=(Slurp const&) -> Slurp& = delete;
};

Being able to very cheaply moved the slurped data around could be very handy indeed.

    [[gnu::nonnull]] Slurp(const char path[]);

Writing functions that look like they take arrays like this is a really bad idea. (Except for main(), which is just special in a lot of ways.) I’ve lost count of the number of newbies who think that this actually takes an array argument, and then think they can do sizeof(path) / sizeof(path[0]) to get the number of elements. Yes, you and I and other experts know what’s really going on, but you shouldn’t be writing code that’s only readable by experts. Down that path lies madness.

On the other hand Slurp(const char* path) or Slurp(char const* path) are unlikely to confuse beginners, and, really, tell the truth of what’s going on much better.

As a general rule, all constructors that take a single argument should be marked explicit (unless you have a damn good reason not to). Yes, explicit should really be the default, but that ship sailed a long time ago.

In addition to this constructor, it would help usability if you could use strings or—even better—actual paths with this class. It would be trivial to do, too:

explicit Slurp(std::string const& path) :
Slurp{path.c_str()}
{}

explicit Slurp(std::filesystem::path const& path) :
Slurp{path.c_str()}
{}

Now this next bit is going to get really complicated. Brace yourself.

### gnu::const on non-static member functions

    [[gnu::const, gnu::returns_nonnull]] const char *begin() const noexcept;

The issue here the gnu::const attribute. I am not a GCC expert, but from my understanding of gnu::const, it sure seems like using it for a member function is an absolutely terrible, and dangerous, idea.

Here’s why.

gnu::const means that a function’s output depends ONLY on the function’s arguments. And I mean ONLY. The function’s behaviour cannot be affected by anything else; it doesn’t matter what state the rest of the program happens to be in (well, I mean, the state has to be valid… if the abstract machine’s state is kerfuckt, then all bets are off in any case), the function’s return value will ALWAYS be the same, given the same arguments.

This makes sense for, say, a function like: auto twice(int i) { return i * 2; }. It doesn’t matter what’s going on with the rest of the program, like how many exceptions are in flight or how many files are open or how much memory has been allocated. So long as the function can complete without UB (for example, there is no stack overflow), then it will ALWAYS return the same result for the same arguments. twice(32) will ALWAYS return 64. ALWAYS. GCC can use that information to optimize. For example, if you call twice(14) a million times in your program, GCC might cache the result of the first call, and then replace every other call with just loading that cached value.

So far none of this sounds problematic, I suppose. Well, now let’s talk about non-static member functions.

There is a very important key feature of non-static class member functions: the implicit this.

You see, every non-static class member function has a hidden this argument. We usually pretend it’s the first argument (but it doesn’t actually have to be, in reality, it’s usually passed in a specific register), so your begin() function actually looks like begin(Slurp const*):

// what you write:
auto Slurp::begin() const noexcept -> char const*
{
return begin_;
}

// what the compiler actually implements:
auto Slurp__begin(Slurp const* this) noexcept -> char const*
{
return this->begin_;
}

Normally, this isn’t something you need to think about.

However… remember what gnu::const means. It means that begin() will ALWAYS return the same result, when called with the same argument.

So let’s imagine you write this:

auto const s = Slurp{"data.file"};

auto p1 = s.begin();
auto p2 = s.begin();
auto p3 = s.begin();
auto p4 = s.begin();
auto p5 = s.begin();

The above is theoretically:

auto const s = Slurp{"data.file"};

auto p1 = Slurp__begin(&s);
auto p2 = Slurp__begin(&s);
auto p3 = Slurp__begin(&s);
auto p4 = Slurp__begin(&s);
auto p5 = Slurp__begin(&s);

But since begin() is gnu::const, once we know the result of calling it with address of s, we can assume that will ALWAYS be the result of calling it with that address. So the compiler could do:

auto p1 = Slurp__begin(&s);
auto p2 = p1;
auto p3 = p1;
auto p4 = p1;
auto p5 = p1;

(And it probably will.)

So far no problems. There is currently no way for the value of begin_ to change once a Slurp object has been constructed, so the above code is a perfectly valid optimization.

However, suppose you implemented moving—or copying, or any way to change the data in a Slurp object. Then you wrote this:

auto const s = Slurp{"data.file"};

auto p1 = s.begin();
auto p2 = s.begin();

s = Slurp{"other.file"};    // this line will construct a new Slurp,
// and REPLACE the internals of s with the
// new data

auto p3 = s.begin();
auto p4 = s.begin();

Now the address of s hasn’t changed… but it’s internal data—in particular, the value of begin_has changed. The compiler sees the above code as:

auto const s = Slurp{"data.file"};

auto p1 = Slurp__begin(&s);
auto p2 = Slurp__begin(&s);

s = Slurp{"other.file"};

auto p3 = Slurp__begin(&s);
auto p4 = Slurp__begin(&s);

And it does the optimization:

auto const s = Slurp{"data.file"};

auto p1 = Slurp__begin(&s);
auto p2 = p1;

s = Slurp{"other.file"};

auto p3 = p1;
auto p4 = p1;

Except after the move assignment, p1 is now a dangling pointer. Which means p3 and p4 are dangling pointers, too.

This is a disaster.

Now you may be thinking: “Meh, no problem; my class is non-copyable, non-movable, and immutable. This problem doesn’t affect me.”

Wanna bet?

Remember, gnu::const isn’t about the object. It’s about the arguments to the function. In the example above, the object’s internals changed, but the address of the object didn’t, so the argument to begin() didn’t either, and that was the problem. So if the object’s internals can’t possibly change… no problem, right?

Wrong.

Because even if it is impossible to change the internals of an object, it is still possible to create an object… destroy it… then create a whole new object at the same memory address… and it is the address that matters. Because that is what is being used as the function argument to begin().

Here’s just one example of how that could happen:

auto s = std::optional<Slurp>{std::in_place, "data.file"};

auto p1 = s->begin();
auto p2 = s->begin();

s.emplace("other.file");    // destroys contained Slurp object, then
// constructs a new one in the same place... so it
// has the same address

auto p3 = s->begin();
auto p4 = s->begin();

Which eventually optimizes to:

auto s = std::optional<Slurp>{std::in_place, "data.file"};

auto p1 = s->begin();
auto p2 = p1;

s.emplace("other.file");

auto p3 = p1;
auto p4 = p1;

And ends up with the same problem.

But you don’t even need std::optional! This problem could be triggered with even this code:

auto f()
{
auto s = Slurp{"data.file"};

auto p = s.begin();
}

auto some_other_function()
{
f();
f();
}

What might happen above is the first call to f() happens to place s at some address on the stack… then when f() ends, s is destroyed… then f() is called again and just happens to place the new s at the same address on the stack. But of course, the allocation done within the Slurp constructor (the posix_memalign() call) may be entirely different. So p SHOULD have two entirely different values with the two successive calls to f() even though the address of s (and thus, the argument passed to begin()) is the same for both calls. But gnu::const says it won’t; gnu::const says that if the argument is the same, the result is the same.

Once again: catastrophe.

To get a clearer view of why, imagine the calls to f() are inlined. Inlining isn’t necessary to trigger the problem, but it helps illustrate. If the calls are inlined, some_other_function() will look like:

auto some_other_function()
{
{
auto s = Slurp{"data.file"};

auto p = s.begin();

// s is destroyed here
}

{
auto s = Slurp{"data.file"};

auto p = s.begin();

// s is destroyed here
}
}

After the first s is destroyed, the second s may very well be constructed at the exact same location… let’s call it STACK_ADDRESS. So:

auto some_other_function()
{
{
auto s = Slurp{"data.file"};

auto p = Slurp__begin(STACK_ADDRESS);

// s is destroyed here
}

{
auto s = Slurp{"data.file"};

auto p = Slurp__begin(STACK_ADDRESS);

// s is destroyed here
}
}

The two calls to begin() are identical, and gnu::const promises that that means the result will be identical, so the compiler might do:

auto some_other_function()
{
auto __p_inlined = static_cast<char const*>(nullptr);

{
auto s = Slurp{"data.file"};

__p_inlined = Slurp__begin(STACK_ADDRESS);
auto p = __p_inlined;

// s is destroyed here
}

{
auto s = Slurp{"data.file"};

auto p = __p_inlined;

// s is destroyed here
}
}

And, boom.

Whew, that was a long digression. The tl;dr is “don’t use gnu::const for (non-static) member functions”… but I think it’s important to understand why.

And, really, it seems kinda pointless in modern C++. If it is possible for a function to be gnu::const, then it is almost certainly also possible for it to be constexpr… in which case the compiler can optimize it far more aggressively than it could possibly optimize a function that is merely gnu::const (but not wholly visible).

### Slurp.cc

using namespace std;

Never, ever do this at file scope. No, not even in a non-header file. Just use std::.

#define ERROR_IF(CONDITION) do {\
if (__builtin_expect((CONDITION), false))\
throw system_error(errno, system_category());\
} while (false)

#define ERROR_IF_POSIX(POSIX_PREFIXED_FUNCTION_CALL) do {\
const int error_number = posix_##POSIX_PREFIXED_FUNCTION_CALL;\
if (__builtin_expect(error_number, 0))\
throw system_error(error_number, system_category());\
} while (false)

Avoid macros. There’s no good reason for them here. Inlined functions will be just as efficient… not that it really matters, because these are error cases after all. For example:

namespace {

auto error_if(int result)
{
if (result != 0) [[unlikely]]
throw std::system_error{errno, std::system_category()};
}

auto error_if_posix(int error_number)
{
if (error_number != 0) [[unlikely]]
throw std::system_error{error_number, std::system_category()};
}

} // anonymous namespace

Slurp::Slurp(char const* path)
{
// ... [snip] ...

error_if_posix(::posix_fadvise(fd, 0, 0, POSIX_FADV_SEQUENTIAL));

struct ::stat file_statistics{};
error_if(::fstat(fd, &file_statistics));

// ... etc.

But even that’s not really great, because it just obfuscates the function’s logic. I don’t see how it’s really an improvement over:

Slurp::Slurp(char const* path)
{
// ... [snip] ...

if (auto const res = ::posix_fadvise(fd, 0, 0, POSIX_FADV_SEQUENTIAL); res != 0) [[unlikely]]
throw std::system_error{res, std::system_category()};

struct ::stat file_statistics{};
if (auto const res = ::fstat(fd, &file_statistics); res != 0) [[unlikely]]
throw std::system_error{errno, std::system_category()};

// ... etc.

The only thing I might do is maybe avoid repeating std::system_category by doing:

namespace {

auto posix_error(int err)
{
return std::system_error{err, std::system_category()};
}

} // anonymous namespace

Slurp::Slurp(char const* path)
{
// ... [snip] ...

if (auto const res = ::posix_fadvise(fd, 0, 0, POSIX_FADV_SEQUENTIAL); res != 0) [[unlikely]]
throw posix_error(res);

auto file_statistics = ::stat{};
if (auto const res = ::fstat(fd, &file_statistics); res != 0) [[unlikely]]
throw posix_error(errno);

// ... etc.

I just don’t see the benefit of hiding the throws. Yes, I get that it is somewhat ugly to have to check the return code every function call… but this is a C API; that’s just the nature of the game.

If you really don’t want all that ugly C-ish code, the right move is to properly wrap the API… not to obfuscate your logic by hiding it behind macros. For example:

Slurp::Slurp(char const* path)
{
auto const fd = posix::open(path, O_RDONLY);

posix::fadvise(fd, 0, 0, POSIX_FADV_SEQUENTIAL);

auto const file_statistics = posix::fstat(fd);

// ... etc.

(Note that implied is that each call above throws on failure, like a good C++ API should.)

But it’s really an all-or-nothing deal. Either completely wrap the API into a proper C++ API, so the usage just looks like regular C++… which maintainers of the code should be familiar with. Or use the C API as-is—don’t hide it behind macros—so people who actually know the POSIX API inside and out can look at your code and immediately spot what’s going on (and any potential problems). Going half-way just annoys both parties, with no real benefits.

    class FileDescriptor {
const int fd_;
public:
FileDescriptor(int fd) : fd_{fd} {ERROR_IF(fd_ == -1);}
~FileDescriptor() {close(fd_);}
[[gnu::const]] operator int() const noexcept {return fd_;}
};

This is a really useful class! Rather than burying it as an implementation detail, it could be part of a public API. You’d need to flesh it out a bit (and make it non-copyable, but perhaps movable), but it’s pretty cool.

FileDescriptor fd{open(path, O_RDONLY)};

This goes back to the “never use using namespace std;” rule: function names like open() are extremely common, and there is a very high chance of conflict between the POSIX function, and some function you might put in your class in the future. If that were to happen, if you are lucky, you’ll get a compile error… but I’ve seen cases where behaviour subtly changes because the wrong function is called.

To make sure you get the right function, you should explicitly qualify it. In this case, because it’s a POSIX function that doesn’t exist in a namespace, you have to use ::open().

You should do this for all the POSIX functions.

constexpr char32_t NullChar = U'\0';

I’m a little baffled by the logic here.

Okay, so you want the data to be null-terminated. Fine, sure. But… you’re treating the data as an array of char… not char32_t. So… huh?

The only thing I can figure is that you’re imagining that maybe someone might possibly want to reinterpret_cast the char const* to char32_t const*, and treat the data as a UTF-32 NUL-terminated string. That seems a pretty odd thing to presume. And it’s even odder when you consider that the Slurp class doesn’t even model a NUL-terminated string in any sense. There is no c_str() function, for example. Even if there were, it’s hard to see why it should return a char32_t const* and not a char const*.

Aside from being odd, it may also be wrong. You are aligning the memory according to file_statistics.st_blksize… not alignof(char32_t). In other words, if you’re thinking that you can just reinterpret_cast the result of begin() to a char32_t const* and treat it as a UTF-32 string… you may be in for a nasty surprise when your program crashes due to a a misaligned read. Or, possibly, just garbled text (that is probably illegal Unicode). (And all this is assuming the endianness of the UTF-32 data is correct.)

I mean, if you really wanted to read file data as a UTF-32 string, why wouldn’t you read it into a std::u32string? Similarly, if you want a regular string, why not read it directly into a std::string? In either case, all the allocation and alignment headaches are taken care of, and the NUL termination is automatic. And if you want binary data, why not read directly into a std::vector<std::byte>?

I’ll be getting back to that idea later.

ERROR_IF_POSIX(memalign(&begin_, blksize, bufsize));

Honestly, I don’t see the point of this.

Okay, I get that the file you’re reading probably has a preferred block size, and if you were mem-mapping the file, or reading it directly from storage, then yes, you’d need aligned memory. But you’re doing neither.

Internally, when you open the file (and specify the access characteristics via posix_fadvise()) a buffer gets created. (I think you can avoid this by using direct input, but that’s beyond the scope here. mmap() is another way to avoid the buffer, but that’s another beast altogether.) Then you are reading from that buffer into your program’s memory space. Now, the buffer needs to be aligned to the block size of the file, so that the file can be read efficiently a block at a time. But your memory… the memory you allocate to read the file data (from the buffer) into… that does not need to be aligned to the block size of the file… because with that memory, you are not reading directly from the file, you are reading memory-to-memory (from the kernel’s buffer, to your array).

Follow? Block-alignment only matters for the initial read of the file from the device. THAT needs to be block-aligned. But unless you’re doing direct input, that’s happening in the kernel (usually), and then read() copies from that buffer to your desired location. Once the data is already in RAM, it’s just a memory-to-memory copy… which is just a rep stosb or some SSE or AVX operation. That doesn’t need to be aligned.⁕

⁕ (Ugh, in order to make this not too absurdly long, I have to gloss over or outright ignore a lot of hyper-technical details, so I know there are nerds out there just champing at the bit to “akshually…” me in the comments. Okay, yes, the data does need to be aligned for super-fast memory-to-memory transfer… both the source and destination data, with penalties for not being aligned varying depending on your processor architecture. HOWEVER… the default C++ allocator in any serious-grade standard library is almost certainly already aligning every allocation to take advantage of those efficiencies by default. So unless you’re doing something “clever”, you don’t need to worry about it.)

In other words, you’re wasting your time using posix_memalign(). You gain nothing from it. You could just as well be using malloc() or new.

Which brings us back again to the idea of reading directly into std::string or std::vector<std::byte>. But let me postpone that for a bit more, because we’re almost at the end of the code, so let’s get through that first.

ERROR_IF_POSIX(madvise(begin_, bufsize, POSIX_MADV_SEQUENTIAL));

I’ll admit that I don’t understand how this call isn’t failing. I thought posix_madvise() was only for memory-mapped memory. I would have guessed you’d get a EBADF from trying to use it on random memory. TIL, I guess.

But I do suspect the call is just being ignored. My understanding of POSIX_MADV_SEQUENTIAL is that it signals to the kernel to make a larger look-ahead buffer… which makes sense if you’re mem-mapping a file into memory, but is just nonsense for RAM. (Where you gonna put the look-ahead buffer for RAM? Super-duper-RAM?)

ERROR_IF(read(fd, begin_, fsize) == -1);

read() can fail spuriously for a number of reasons, like a signal interrupting it, and end up doing only a partial read. You can’t just write a read() call like this; you need to call it in a loop and keep doing reads until you’ve read fsize bytes in total or until read() returns zero.

memcpy(end_, &NullChar, sizeof NullChar); // ensure null termination

I mean… why not just std::fill_n(end, 4, 0);? It would probably be faster. (Unless your optimizer is smart enough to transform that memcpy() into a memset(), in which case it will probably be equally fast.)

## Reading directly into containers

Okay, the biggest problem with Slurp is that it’s a specialized type. Even if it had a perfect interface, it’s still not a string. If I need a std::string, which is very common, then I’ll be forced to copy the entire Slurp object into a std::string. What a waste.

What if you just read the file’s contents directly into a std::string. Or, if it’s binary data, into a std::vector<std::byte>. Remember, there is no need to use posix_memalign() to specially align the memory to the original file’s preferred block size. The default allocator will align things just fine. So you can just use standard containers.

(But, FYI, even if the default allocator didn’t align things just fine, it would still be a better solution to use a custom allocator, and standard containers, rather that rolling your own, half-assed container.)

For example, to read into a string (note: untested, rough code):

auto read_into_string(std::filesystem::path const& path)
{
// file_descriptor is a class basically identical to yours
auto const fd = file_descriptor{::open(path.c_str(), O_RDONLY)};

// get the file size from the file descriptor
auto const file_size = [&]
{
struct ::stat fstats{};
if (auto const res = ::fstat(fd, &fstats); res != 0)
throw std::system_error{errno, std::system_category()};
return fstats.st_size;
}();

// create the string we'll be reading into
auto s = std::string{};
s.resize(file_size);

// this is the read loop
auto pos = std::string::size_type{};
while (true)
{
auto const last_read = ::read(fd, s.data() + pos, s.size() - pos);
if (last_read < 0)
throw std::system_error{errno, std::system_category()};
else if (last_read == 0)
break;

// this should never happen
if (last_read > (s.size() - pos)) [[unlikely]]
throw std::runtime_error{"more data in file than file size"};

pos += last_read;
}

// this should never happen
if (pos != s.size()) [[unlikely]]
throw std::runtime_error{"failed to read entire file"};

return s;
}

And the usage:

auto const string_data = read_into_string("file.txt");

And string_data is a std::string.

The only thing that I suspect might make this code slower than your code, and only for very, very large files, is that there is no way (AFAIK) to allocate a std::string without initializing its contents. In other words, s.resize(n) doesn’t does just allocate n bytes (internally), it also zeros those bytes… which is a total waste of time for us. You shouldn’t have this inefficiency if you’re reading binary data into a std::vector<std::byte> (for which the code is otherwise identical).

But even that std::string inefficiency should be tolerable (unless your files are huge), because you get the benefit of actually having an honest-to-goodness string. Which basically works everywhere, and which many APIs specifically need, in one form or another.

And, frankly, I seriously doubt you’ll even notice it, because, again, unless you’re dealing with huge files, the cost of that zeroing will be completely dwarfed by the I/O costs.

Anywho, this is getting way too long, and I’ve been sitting on it for days, so, I’ll just go ahead and publish what I’ve got.

• Thanks for giving me such a thorough answer! My detailed response to you is in my own answer below. May 13 at 1:05

I would expect Slurp.cc to include Slurp.h - in fact, I can't see how you could compile it without the definitions there. I recommend making it the first include, to detect whether any dependencies are introduced.

using namespace std;

Bringing all of the std identifiers into the global namespace is usually a bad idea - just write std::size_t and other identifiers in full.

It seems strange to use a std::char32_t as terminator for a string of char. If we need to read files with different character types, then I would expect to handle that using a template class.

Obviously files with embedded null characters will appear truncated when interpreted as (null-terminated) C strings; it's good that you provide begin() and end() so users can wrap in std::string_view - although start and length would feed better to its constructor.

• 1) In case it wasn't clear, my code snippet was supposed to be Slurp.cc after Slurp.h was included. 2) It's just char32_t, no std:: (it's a fundamental type). 3) The character encoding (or lack thereof) of a file is not necessarily known at compile-time so I find it most correct to always handle the case where the file is encoded as UTF-32. In "real" code Slurp would provide different "views" or iterators for different character encodings. May 12 at 20:58

Here is how I've updated my code, for C++20, to make the API more usable and full-featured, and with feedback from indi. After that I've written a response to indi's feedback, since I think he deserves one and it's too long and detailed to do so in comments.

Demangle.h

#ifndef DEMANGLE_H
#define DEMANGLE_H

#include <cxxabi.h>

#include <cstdlib>
#include <new>
#include <string_view>
#include <typeinfo>

template<typename T>
class Demangle {
char *name{nullptr};

Demangle(std::string_view s)
: name{static_cast<char*>(std::calloc(s.size() + 1, 1))}
{
if (!name)
throw std::bad_alloc{};

s.copy(name, s.size());
}
public:
Demangle() {
int status;
name = abi::__cxa_demangle(
typeid(T).name(), nullptr, nullptr, &status);

if (!name)
throw std::bad_alloc{};
}

Demangle(const Demangle &other) : Demangle{std::string_view{other}} {}

Demangle(Demangle &&other) noexcept {std::swap(name, other.name);}

Demangle &operator=(Demangle rhs) noexcept {
std::swap(name, rhs.name);
return *this;
}

~Demangle() {std::free(name);}

operator std::string_view() const noexcept {return name;}
};

#endif//DEMANGLE_H

FileDescriptor.h

#ifndef FILE_DESCRIPTOR_H
#define FILE_DESCRIPTOR_H

class FileDescriptor {
const int fd_;
public:
explicit FileDescriptor(int fd);
[[gnu::nonnull]] FileDescriptor(const char path[], int oflag);
/* std::string / std::filesystem::path / similar constructor
*/
template<class T>
FileDescriptor(const T &t, int oflag)
: FileDescriptor{t.c_str(), oflag} {}

~FileDescriptor();

[[gnu::const]] operator const int&() const noexcept;
};

#endif//FILE_DESCRIPTOR_H

ProjectConcepts.h

#ifndef PROJECT_CONCEPTS_H
#define PROJECT_CONCEPTS_H

#include <type_traits>

namespace Project {

template<typename T>
concept TriviallyCopyable = std::is_trivially_copyable_v<T>;

}//Project

#endif//PROJECT_CONCEPTS_H

Slurp.h

#ifndef SLURP_H
#define SLURP_H

#include "Demangle.h"
#include "FileDescriptor.h"
#include "ProjectConcepts.h"

#include <filesystem>
#include <initializer_list>
#include <span>
#include <stdexcept>
#include <string>

class Slurp {
template<typename T>
static constexpr bool IsSmol = sizeof(T) == 1;

struct Bounds {void *begin{nullptr}, *end;} b;

template<Project::TriviallyCopyable T>
std::initializer_list<T*> init() const noexcept(IsSmol<T>) {
if (size() % sizeof(T)) [[unlikely]] {
const Demangle<T> tp_name;
throw std::logic_error{std::string{
"can't interpret file as array of "} + tp_name};
}
return {static_cast<T*>(b.begin), static_cast<T*>(b.end)};
}
public:
Slurp(const Slurp&) = delete;
Slurp &operator=(const Slurp&) = delete;

Slurp(Slurp &&other) noexcept;
Slurp &operator=(Slurp &&rhs) noexcept;

explicit Slurp(const FileDescriptor &fd);
[[gnu::nonnull]] Slurp(const char path[]);
explicit Slurp(const std::string &path);
explicit Slurp(const std::filesystem::path &path);

~Slurp();

[[gnu::returns_nonnull]] const char *begin() const noexcept;
[[gnu::returns_nonnull]] const char *end() const noexcept;

std::size_t size() const noexcept;

template<Project::TriviallyCopyable T, class Traits=std::char_traits<T>>
operator std::basic_string_view<T, Traits>() const noexcept(IsSmol<T>) {
return init<T>();
}
template<Project::TriviallyCopyable T>
operator std::span<T>() & noexcept(IsSmol<T>) {
return init<T>();
}
template<Project::TriviallyCopyable T>
operator std::span<const T>() const noexcept(IsSmol<T>) {
return init<T>();
}
};

#endif//SLURP_H

SystemErrorIf.h

#ifndef SYSTEM_ERROR_IF_H
#define SYSTEM_ERROR_IF_H

#include <cerrno>
#include <system_error>

namespace SystemErrorIf {

constexpr void condition(bool condition) {
if (condition) [[unlikely]]
throw std::system_error{errno, std::system_category()};
}

constexpr void code(int error_number) {
if (error_number) [[unlikely]]
throw std::system_error{error_number, std::system_category()};
}

}//SystemErrorIf

#endif//SYSTEM_ERROR_IF_H

FileDescriptor.cc

#include "FileDescriptor.h"
#include "SystemErrorIf.h"

#include <fcntl.h>
#include <unistd.h>

FileDescriptor::FileDescriptor(int fd) : fd_{fd} {
SystemErrorIf::condition(fd_ == -1);
}
FileDescriptor::FileDescriptor(const char path[], int oflag)
: FileDescriptor{::open(path, oflag)} {}

FileDescriptor::~FileDescriptor() {::close(fd_);}

FileDescriptor::operator const int&() const noexcept {return fd_;}

Slurp.cc

#include "Slurp.h"
#include "SystemErrorIf.h"

#include <fcntl.h>
#include <sys/mman.h>
#include <sys/stat.h>
#include <unistd.h>

#include <cstring>

Slurp::Slurp(Slurp &&other) noexcept {std::swap(b, other.b);}

Slurp &Slurp::operator=(Slurp &&rhs) noexcept {
std::swap(b, rhs.b);
return *this;
}

Slurp::Slurp(const FileDescriptor &fd) {
namespace ErrorIf = SystemErrorIf;

ErrorIf::code(::posix_fadvise(fd, 0, 0, POSIX_FADV_SEQUENTIAL));

struct ::stat file_statistics;

ErrorIf::condition(::fstat(fd, &file_statistics));

::ssize_t fsize = file_statistics.st_size;

const auto bufsize = static_cast<std::size_t>(fsize) + sizeof(char32_t),
blksize = static_cast<std::size_t>(file_statistics.st_blksize);

ErrorIf::code(::posix_memalign(&b.begin, blksize, bufsize));

ErrorIf::code(::posix_madvise(b.begin, bufsize, POSIX_MADV_SEQUENTIAL));

b.end = b.begin;

while (fsize > 0) {
const auto nread = ::read(
fd, b.end, static_cast<std::size_t>(fsize));

ErrorIf::condition(nread == -1);
b.end = static_cast<char*>(b.end) + nread;
fsize -= nread;
}

std::memset(b.end, 0, sizeof(char32_t)); // ensure null termination
}
Slurp::Slurp(const char path[]) : Slurp{FileDescriptor{path, O_RDONLY}} {}
Slurp::Slurp(const std::string &path) : Slurp{path.c_str()} {}
Slurp::Slurp(const std::filesystem::path &path) : Slurp{path.c_str()} {}

Slurp::~Slurp() {std::free(b.begin);}

const char *Slurp::begin() const noexcept {return static_cast<char*>(b.begin);}
const char *Slurp::end() const noexcept {return static_cast<char*>(b.end);}

std::size_t Slurp::size() const noexcept {return end() - begin();}

indi, first of all, thank you for putting so much effort into your answer. Even if I don't agree with everything you said, I really appreciate the work you put in.

The convention in C++ is to put all the type information together

I don't think there's any particular "convention" in C++ (or C for that matter). It's my personal preference/style to put sigils (*, &, &&) next to identifiers and not types.

Rather than holding two pointers, it might make more sense to hold a pointer and a size.

Again, this is more of a personal preference thing.

I also don’t get why [the begin_ and end_ pointers] are void*

Two reasons, first because the file may have any or no character encoding, so void* is semantically appropriate (whereas char* implies ASCII-compatible text and unsigned char* or std::byte* imply an unstructured stream of binary data), and second because posix_memalign(3) (which, as you say, may be overkill) takes a void**, to which a char** would need to be reinterpret_cast, and I don't want to open up that whole can of worms.

This [deleted default constructor] is technically unnecessary

Good point.

if your intention is to disable copying, then [the deleted copy constructor and copy assignment operator] are malformed.

I don't believe so, but it's not worth looking into this just to save a few characters, so I'll write these your way from now on.

you might consider making the type movable.

Done.

Writing functions that look like they take arrays like this is a really bad idea.

This is a low-level I/O routine, it's not written for beginners. While it does go back to a really bad language design decision on the part of the creators of C, I don't think you have to be an "expert" to know that foo(const char bar[]) takes a pointer. I think it's more semantic to write *bar when bar points to a single object and bar[] when, as in this case, it points to an array.

I am flattered you think I'm an expert! :)

As a general rule, all constructors that take a single argument should be marked explicit (unless you have a damn good reason not to)

I like explicit too but in this case I don't think there's any need for extra ceremony if a class has an operator char*/const char*.

In addition to this constructor, it would help usability if you could use strings or—even better—actual paths with this class.

Done, and part of what I intended to do anyway.

gnu::const on non-static member functions

I don't think your interpretation of gnu::const is correct. The compiler knows when a new object has been created, even if it shares the same address as a previous object. But it's a moot point since as per your suggestion I made Slurp movable, and so removed gnu::const from its member functions.

Never, ever [using namespace std;] at file scope. No, not even in a non-header file.

Opinions on this differ, I don't think it's that big a deal in a non-header file, particularly a short one with few includes, but I'll remove it anyway.

Avoid macros.

I was targeting C++11 in the original question's code, and wanted to ensure these were inlined. Since I'm targeting C++20 in this code I can and have made them constexpr functions instead.

But even that’s not really great, because it just obfuscates the function’s logic. I don’t see how it’s really an improvement

I guess it's a bit of a Ruby-ism (Don't Repeat Yourself) that I tend to do in other languages also. When I see repeated logic I extract it to a function (or macro as the case may be).

If you really don’t want all that ugly C-ish code, the right move is to properly wrap the API…

Maybe, but I'd say that's out of scope for now.

This is a really useful class! Rather than burying it as an implementation detail, it could be part of a public API.

Done.

You’d need to flesh it out a bit (and make it non-copyable, but perhaps movable)

I don't think it makes sense for a file descriptor to be movable; a POSIX file descriptor is conceptually unique, immutable and identifies a single object which is opened and closed exactly once (per file descriptor).

To make sure you get the right function, you should explicitly qualify it.

Done, although IIRC the POSIX standard may allow for some of these functions to be implemented as macros, in which case this may lead to compilation errors on certain (hypothetical) platforms.

The only thing I can figure is that you’re imagining that maybe someone might possibly want to reinterpret_cast the char const* to char32_t const*, and treat the data as a UTF-32 NUL-terminated string.

Yes, that's what I'm intending.

And it’s even odder when you consider that the Slurp class doesn’t even model a NUL-terminated string in any sense.

It was intentionally bare-bones for explication purposes, as you can see it's more full-featured in my updated code.

You are aligning the memory according to file_statistics.st_blksize… not alignof(char32_t). In other words, if you’re thinking that you can just reinterpret_cast the result of begin() to a char32_t const* and treat it as a UTF-32 string… you may be in for a nasty surprise when your program crashes due to a a misaligned read.

Well I'm assuming (perhaps strictly incorrectly) that st_blksize (512 on Linux) > alignof(char32_t) (4), and I'm pretty sure any allocation function (posix_memalign(3) included) will always allocate memory that's at least suitably aligned for std::max_align_t.

And all this is assuming the endianness of the UTF-32 data is correct.

Probably something I should look into, though I haven't addressed it in the below code.

if you want a regular string, why not read it directly into a std::string?

Because it introduces overhead and I intend for this class to be as efficient as possible. It's trivial to construct a std::string (or std::u32string, std::vector<std::byte>, etc.) from it if desired.

Honestly, I don’t see the point of [using posix_memalign(3)].

You may be correct, I need to write some benchmarks to see what the performance of posix_memalign(3) is compared to malloc(3).

I thought posix_madvise() was only for memory-mapped memory.

I don't think so? The man page says nothing about that. Again I should probably write a benchmark.

read() can fail spuriously for a number of reasons, like a signal interrupting it, and end up doing only a partial read. You can’t just write a read() call like this; you need to call it in a loop and keep doing reads until you’ve read fsize bytes in total

Good point, and done.

I mean… why not just std::fill_n(end, 4, 0);? It would probably be faster. (Unless your optimizer is smart enough to transform that memcpy() into a memset(), in which case it will probably be equally fast.)

Good point re: memset(3). In this ("C with classes"-ish) context I prefer it to std::fill_n().

Reading directly into containers

I largely agree with you, but I'm not doing so for three reasons: 1) Part of the point of this exercise was to be as efficient as possible, and reading into a container creates overhead. 2) Slurp is supposed to model any possible character encoding (or lack thereof) that a file might have. 3) C++17 and C++20 eliminate much of the need for this (at least for newly-written APIs) with std::basic_string_view and std::span.


text:
Numerical Integration of Newton's Equation of Motion

posted January 30, 2021
by H. Gould, J. Tobochnik, and W. Christian

Appendix 3B of An Introduction to Computer Simulation Methods 3rd Ed., by H. Gould, J. Tobochnik, and W. Christian (2007).

Download all of Chapter 3 as PDF from the Open Source Physics collection

View/Download other chapters from this book from the Open Source Physics collection

We summarize several of the common finite difference methods for the solution of Newton's equations of motion with continuous force functions. The number and variety of algorithms currently in use is evidence that no single method is superior under all conditions.

To simplify the notation, we consider the motion of a particle in one dimension and write Newton's equations of motion in the form \label{eq:motion/eqnewton} \begin{align} {dv \over dt} &= a(t),\\ {dx \over dt} &= v(t), \end{align} where $a(t) \equiv a\left(x(t),v(t),t\right)$. The goal of finite difference methods is to determine the values of $x_{n+1}$ and $v_{n+1}$ at time $t_{n+1} = t_n + dt$. We already have seen that $dt$ must be chosen so that the integration method generates a stable solution. If the system is conservative, $dt$ must be sufficiently small so that the total energy is conserved to the desired accuracy.

The nature of many of the integration algorithms can be understood by expanding $v_{n+1} = v(t_n + dt)$ and $x_{n+1} = x(t_n + dt)$ in a Taylor series. We write \begin{align} v_{n+1} &= v_n + a_n dt + O \big ((dt)^2\big), \label{eq:motion/taylorx} \\ x_{n+1} &= x_n + v_n dt + \frac{1}{2} a_n (dt)^2 + O \left((dt)^3\right). \label{eq:motion/taylorv} \end{align} The familiar Euler algorithm is equivalent to retaining the $O(dt)$ terms in Eqs. (\ref{eq:motion/taylorx}) and (\ref{eq:motion/taylorv}): \label{eq:motion/firstorder} \begin{align} v_{n+1} &= v_n + a_n dt \\ x_{n+1} &= x_n + v_n dt. \qquad \mbox{(Euler algorithm)} \end{align} Because order $dt$ terms are retained in (\ref{eq:motion/firstorder}), the local truncation error, the error in one time step, is order $(dt)^2$. The global error, the total error over the time of interest, due to the accumulation of errors from step to step is order $dt$. This estimate of the global error follows from the fact that the number of steps into which the total time is divided is proportional to $1/dt$. Hence, the order of the global error is reduced by a factor of $1/dt$ relative to the local error. We say that an algorithm is $n$th order if its global error is order $(dt)^n$. The Euler algorithm is an example of a *first-order* algorithm.

The Euler algorithm is asymmetrical because it advances the solution by a time step $dt$, but uses information about the derivative only at the beginning of the interval. We already have found that the accuracy of the Euler algorithm is limited and that frequently its solutions are not stable. We also found that a simple modification of (\ref{eq:motion/firstorder}) yields solutions that are stable for oscillatory systems. For completeness, we repeat the Euler-Cromer algorithm here: \begin{align} v_{n+1} &= v_n + a_n dt, \\ x_{n+1} &= x_n + v_{n+1}dt. \qquad \mbox{(Euler-Cromer algorithm)} \end{align}

An obvious way to improve the Euler algorithm is to use the mean velocity during the interval to obtain the new position. The corresponding *midpoint* algorithm can be written as \begin{align} v_{n+1} &= v_n + a_n dt,\quad\mbox{and} \label{eq:motion/midv}\\ x_{n+1} &= x_n + \frac{1}{2} \left(v_{n+1} + v_n\right) dt. \qquad \mbox{(midpoint algorithm)} \label{eq:motion/midx} \end{align} Note that if we substitute (\ref{eq:motion/midv}) for $v_{n+1}$ into (\ref{eq:motion/midx}), we obtain $$x_{n+1} = x_n + v_n dt + \frac{1}{2} a_n \,dt^2.$$ Hence, the midpoint algorithm yields second-order accuracy for the position and first-order accuracy for the velocity. Although the midpoint approximation yields exact results for constant acceleration, it usually does not yield much better results than the Euler algorithm. In fact, both algorithms are equally poor, because the errors increase with each time step.

A higher order algorithm whose error is bounded is the *half-step* algorithm. In this algorithm the average velocity during an interval is taken to be the velocity in the middle of the interval. The half-step algorithm can be written as \begin{align} v_{n + \frac{1}{2}} &= v_{n-\frac{1}{2}} + a_n dt, \label{eq:motion/halfv}\\ x_{n+1} &= x_n + v_{n+\frac{1}{2}} dt. \qquad \qquad \mbox{(half-step algorithm)} \end{align}

Note that the half-step algorithm is not self-starting, that is, (\ref{eq:motion/halfv}) does not allow us to calculate $v_\frac{1}{2}$. This problem can be overcome by adopting the Euler algorithm for the first half step:

$$v_\frac{1}{2} = v_0 + \frac{1}{2} a_0 dt.$$

Because the half-step algorithm is stable, it is a common textbook algorithm. The Euler-Richardson algorithm, a widely used half-step algorithm, can be motivated as follows. We first write $x(t+dt)$ as $$x_1 \approx x(t+dt) = x(t) + v(t)dt + \frac{1}{2} a(t) (dt)^2. \label{eq:motion/taylor11}$$ The notation $x_1$ implies that $x(t+dt)$ is related to $x(t)$ by one time step. We also may divide the step $dt$ into half steps and write the first half step, $x(t+\frac{1}{2} dt)$, as $$\label{eq:motion/taylor2} x(t+\frac{1}{2}dt) \approx x(t) + v(t){dt \over 2}+ \frac{1}{2} a(t) \big ({dt \over 2} \big)^2.$$ The second half step, $x_2(t+dt)$, may be written as $$\label{eq:motion/taylor3} x_2(t+dt) \approx x(t+\frac{1}{2} dt) + v(t+\frac{1}{2} dt){dt \over 2} + \frac{1}{2} a(t + \frac{1}{2} dt)\big ({dt\over 2}\big)^2.$$ We substitute (\ref{eq:motion/taylor2}) into (\ref{eq:motion/taylor3}) and obtain $$x_2(t+dt) \approx x(t) + \frac{1}{2} \big [v(t) + v(t+\frac{1}{2} dt) \big] dt + \frac{1}{2} \big [a(t) + a(t + \frac{1}{2} dt) \big]\big (\frac{1}{2} dt\big)^2. \label{eq:motion/taylor4}$$ Now $a(t+\frac{1}{2} dt) \approx a(t) + \frac{1}{2} a'(t)dt$. Hence to order $(dt)^2$, (\ref{eq:motion/taylor4}) becomes $$x_2(t+dt) = x(t) + \frac{1}{2} \big [v(t) + v(t+\frac{1}{2}dt)\big]dt + \frac{1}{2} \big [2a(t) \big] \big (\frac{1}{2}dt\big)^2. \label{eq:motion/taylor5}$$

We can find an approximation that is accurate to order $(dt)^3$ by combining (\ref{eq:motion/taylor11}) and (\ref{eq:motion/taylor5}) so that the terms to order $(dt)^2$ cancel. The combination that works is $2x_2 - x_1$, which gives the Euler-Richardson result: $$x_{\rm er}(t +dt) \equiv 2x_2(t+dt)- x_1(t+dt) = x(t) + v(t+\frac{1}{2}dt)dt + O(dt)^3.\label{eq:motion/ermy}$$ The same reasoning leads to an approximation for the velocity accurate to $(dt)^3$ giving $$v_{\rm er} \equiv 2v_2(t+dt)- v_1(t+dt) = v(t) + a(t+\frac{1}{2}dt)dt + O(dt)^3. \label{eq:motion/ermv}$$

A bonus of the Euler-Richardson algorithm is that the quantities $|x_2-x_1|$ and $|v_2-v_1|$ give an estimate for the error. We can use these estimates to change the time step so that the error is always within some desired level of precision. We will see that the Euler-Richardson algorithm is equivalent to the second-order Runge-Kutta algorithm (see (\ref{eq:motion/rktwo})).

One of the most common energy drift-free higher order algorithms is commonly attributed to Verlet. We write the Taylor series expansion for $x_{n-1}$ in a form similar to (\ref{eq:motion/taylorx}): $$x_{n-1} = x_n - v_n dt + \frac{1}{2} a_n (dt)^2. \label{eq:motion/forward}$$ If we add the forward and reverse forms, (\ref{eq:motion/taylorx}) and (\ref{eq:motion/forward}) respectively, we obtain $$x_{n+1} + x_{n-1} = 2x_n + a_n (dt)^2 + O \big ((dt)^4 \big)$$ or $$x_{n+1} = 2x_n - x_{n-1} + a_n (dt)^2. \qquad \mbox{(leapfrog algorithm)} \label{eq:motion/verletx}$$

Similarly, the subtraction of the Taylor series for $x_{n+1}$ and $x_{n-1}$ yields

$$v_n = {x_{n+1} - x_{n-1} \over 2dt}. \qquad \qquad \mbox{(leapfrog algorithm)} \label{eq:motion/verletv}$$ Note that the global error associated with the leapfrog algorithm is third-order for the position and second-order for the velocity. However, the velocity plays no part in the integration of the equations of motion. The leapfrog algorithm is also known as the explicit central difference algorithm. Because this form of the leapfrog algorithm is not self-starting, another algorithm must be used to obtain the first few terms. An additional problem is that the new velocity (\ref{eq:motion/verletv}) is found by computing the difference between two quantities of the same order of magnitude. Such an operation results in a loss of numerical precision and may give rise to roundoff errors.

A mathematically equivalent version of the leapfrog algorithm is given by \begin{align} x_{n+1} &= x_n + v_n dt + \frac{1}{2} a_n (dt)^2 \label{eq:motion/velverletx} \\ v_{n+1} &= v_n + \frac{1}{2} \left(a_{n+1} + a_n\right) dt. \qquad \mbox{(velocity Verlet algorithm)} \label{eq:motion/velverletv} \end{align} We see that Eqs. (\ref{eq:motion/velverletx}) and (\ref{eq:motion/velverletv}), known as the *velocity* form of the Verlet algorithm, is self-starting and minimizes roundoff errors. We will refer to Eqs. (\ref{eq:motion/velverletx}) and (\ref{eq:motion/velverletv}) as the Verlet algorithm.

We can derive Eqs. (\ref{eq:motion/velverletx}) and (\ref{eq:motion/velverletv}) from the leapfrog algorithm by the following considerations. We first solve (\ref{eq:motion/verletv}) for $x_{n-1}$ and write $x_{n-1}=x_{n+1} - 2v_n dt$. If we substitute this expression for $x_{n-1}$ into (\ref{eq:motion/verletx}) and solve for $x_{n+1}$, we find the form (\ref{eq:motion/velverletx}). Then we use (\ref{eq:motion/verletv}) to write $v_{n+1}$ as: $$v_{n+1} = {x_{n+2} - x_n \over 2 dt}, \label{eq:motion/vderiv}$$ and use (\ref{eq:motion/verletx}) to obtain $x_{n+2}=2x_{n+1} - x_n + a_{n+1} (dt)^2$. If we substitute this form for $x_{n+2}$ into (\ref{eq:motion/vderiv}), we obtain $$v_{n+1}= {x_{n+1} - x_n \over dt} + \frac{1}{2} a_{n+1} dt. \label{eq:motion/vinterm}$$ Finally, we use (\ref{eq:motion/velverletx}) for $x_{n+1}$ to eliminate $x_{n+1} - x_n$ from (\ref{eq:motion/vinterm}); after some algebra we obtain the desired result (\ref{eq:motion/velverletv}).

Another useful algorithm that avoids the roundoff errors of the leapfrog algorithm is due to Beeman and Schofield. We write the *Beeman* algorithm in the form: \label{eq:motion/beeman} \begin{align} x_{n+1} &= x_n + v_n dt + {1 \over 6}(4a_n -a_{n-1}) (dt)^2 \\ v_{n+1} &= v_n + {1 \over 6} (2a_{n+1} + 5a_n - a_{n-1}) dt. \qquad \mbox{(Beeman algorithm)} \end{align} Note that (\ref{eq:motion/beeman}) does not calculate particle trajectories more accurately than the Verlet algorithm. Its advantage is that it generally does a better job of maintaining energy conservation. However, the Beeman algorithm is not self-starting.

The most common finite difference method for solving ordinary differential equations is the *Runge-Kutta* method. To explain the many algorithms based on this method, we consider the solution of the first-order differential equation $${dx \over dt} = f(x,t). \label{eq:motion/rk}$$ Runge-Kutta algorithms evaluate the rate, $f(x,t)$, multiple times in the interval $[t, t+\Delta t]$. For example, the classic fourth-order Runge-Kutta algorithm, which we will discuss in the following, evaluates $f(x,t)$ at four times $t_n$, $t_n +a_1\Delta t$, $t_n +a_2\Delta t$, and $t_n +a_3 \Delta t$. Each evaluation of $f(x,t)$ produces a slightly different rate $r_1$, $r_2$, $r_3$, and $r_4$. The idea is to advance the solution using a weighted average of the intermediate rates: $$y_{n+1} = y_n + (c_1 r_1 + c_2 r_2 +c_3 r_3 + c_4 r_4)\Delta t.$$

The various Runge-Kutta algorithms correspond to different choices for the constants $a_i$ and $c_i$. These algorithms are classified by the number of intermediate rates $\{r_i\}$, where $i=1, \cdots, N$. The determination of the Runge-Kutta coefficients is difficult for all but the lowest order methods, because the coefficients must be chosen to cancel as many terms in the Taylor series expansion of $f(x,t)$ as possible. The first non-zero expansion coefficient determines the order of the Runge-Kutta algorithm. Fortunately, these coefficients are tabulated in most numerical analysis books.

To illustrate how the various sets of Runge-Kutta constants arise, consider the case $N=2$. The second-order Runge-Kutta solution of (\ref{eq:motion/rk}) can be written using standard notation as: \label{eq:motion/rktwo} \begin{align} x_{n+1} &= x_n + k_2 + O \left( (dt)^3 \right), \quad\mbox{where}\\ k_2 &= f\left(x_n + {k_1 \over 2}, t_n + {dt \over 2}\right) dt. \\ k_1 &= f\left(x_n,t_n\right) dt \end{align} Note that the weighted average uses $c_1=0$ and $c_2=1$. The interpretation of (\ref{eq:motion/rktwo}) is as follows. The Euler algorithm assumes that the slope $f(x_n,t_n)$ at $(x_n,t_n)$ can be used to extrapolate to the next step, that is, $x_{n+1} = x_n + f(x_n,t_n) dt$. A plausible way of making a a more accurate determination of the slope is to use the Euler algorithm to extrapolate to the midpoint of the interval and then to use the midpoint slope across the full width of the interval. Hence, the Runge-Kutta estimate for the rate is $f(x^{\*}, t_n + \frac{1}{2} dt)$, where $x^{\*} = x_n + \frac{1}{2} f(x_n,t_n) dt$ (see (\ref{eq:motion/rktwo})).

The application of the second-order Runge-Kutta algorithm to Newton's equation of motion (\ref{eq:motion/eqnewton}) yields \label{eq:motion/rk2Newton} \begin{align} k_{1v}& = a_n(x_n,v_n,t_n) dt \\ k_{1x}& = v_n dt \\ k_{2v}& = a\left(x_n + {k_{1x} \over 2},v_n + { k_{1v} \over 2}, t + {dt \over 2}\right)dt \\ k_{2x}& = \left(v_n + {k_{1v} \over 2}\right) dt, \end{align} and \label{eq:motion/rk2} \begin{align} v_{n+1}& = v_n + k_{2v} \\ x_{n+1}& = x_n + k_{2x}. \qquad \mbox{(second-order Runge Kutta)} \end{align} Note that the second-order Runge-Kutta algorithm in (\ref{eq:motion/rk2Newton}) and (\ref{eq:motion/rk2}) is identical to the Euler-Richardson algorithm.

Other second-order Runge-Kutta type algorithms exist. For example, if we set $c_1=c_2=\frac{1}{2}$ we obtain the endpoint method: \begin{align} y_{n+1} &= y_n + \frac{1}{2} k_1+ \frac{1}{2} k_2, \quad\mbox{where}\\ k_{1}& = f\left(x_n,t_n\right) dt \\ k_{2}& = f\left(x_n+k_1,t_n+dt\right) dt. \end{align} And if we set $c_1=\frac{1}{3}$ and $c_2=\frac{2}{3}$, we obtain Ralston's method: \begin{align} y_{n+1} &= y_n + \frac{1}{3} k_1+ \frac{2}{3} k_2, \quad\mbox{where}\\ k_{1}& = f\left(x_n,t_n\right) dt \\ k_{2}& = f\left(x_n+\frac{3}{4}k_1,t_n+\frac{3}{4}dt\right) dt. \end{align} Note that Ralston's method does not calculate the rate at uniformly spaced subintervals of $dt$. In general, a Runge-Kutta method adjusts the partition of $dt$ as well as the constants $a_i$ and $c_i$ so as to minimize the error.

In the *fourth-order* Runge-Kutta algorithm, the derivative is computed at the beginning of the time interval, in two different ways at the middle of the interval, and again at the end of the interval. The two estimates of the derivative at the middle of the interval are given twice the weight of the other two estimates. The algorithm for the solution of (\ref{eq:motion/rk}) can be written in standard notation as \label{eq:motion/rkfour} \begin{align} k_1&= f\left(x_n,t_n\right) dt \\ k_2& = f\left(x_n + {k_1 \over 2},t_n + {dt \over 2}\right) dt \\ k_3& = f\left(x_n+{k_2 \over 2},t_n + {dt \over 2}\right) dt \\ k_4& = f\left(x_n + k_3,t_n +dt\right) dt, \end{align} and $$x_{n+1}= x_n + {1 \over 6}(k_1 + 2k_2 + 2k_3 + k_4).$$

The application of the fourth-order Runge-Kutta algorithm to Newton's equation of motion (\ref{eq:motion/eqnewton}) yields \begin{align} k_{1v}&= a\left(x_n,v_n,t_n\right) dt \\ k_{1x}& = v_n dt \\ k_{2v}& = a\left(x_n+{k_{1x} \over 2},v_n+{ k_{1v} \over 2},t_n+{ dt \over 2}\right) dt \\ k_{2x}& = \left(v_n + {k_{1v} \over 2}\right) dt \\ k_{3v}& = a\left(x_n + {k_{2x} \over 2}, v_n + {k_{2v} \over 2}, t_n + {dt \over 2}\right) dt \\ k_{3x}& = \left(v_n + {k_{2v} \over 2}\right) dt \\ k_{4v}& = a\left(x_n + k_{3x}, v_n+ k_{3v}, t+dt\right) dt \\ k_{4x}& = \left(v_n + k_{3v}\right) dt, \end{align} and \begin{align} v_{n+1}& = v_n + {1 \over 6}(k_{1v} + 2k_{2v}+2k_{3v}+k_{4v}) \\ x_{n+1}& = x_n + {1 \over 6}(k_{1x} + 2k_{2x}+2k_{3x}+k_{4x}). \quad \mbox{(fourth-order Runge-Kutta)} \end{align} Because Runge-Kutta algorithms are self-starting, they are frequently used to obtain the first few iterations for an algorithm that is not self-starting.

As we have discussed, one way to determine the accuracy of a solution is to calculate it twice with two different values of the time step. One way to make this comparison is to choose time steps $dt$ and $dt/2$ and compare the solution at the desired time. If the difference is small, the error is assumed to be small. This estimate of the error can be used to adjust the value of the time step. If the error is too large, than the time step can be halved. And if the error is much less than the desired value, the time step can be increased so that the program runs faster.

A better way of controlling the step size was developed by Fehlberg who showed that it is possible to evaluate the rate in such as way as to simultaneously obtain two Runge-Kutta approximations with different orders. For example, it is, possible to run a fourth-order and fifth-order algorithm in tandem by evaluating five rates. We thus obtain different estimates of the true solution using different weighed averages of these rates: \begin{align} y_{n+1} &= y_n + c_1 k_1 + c_2 k_2 +c_3 k_3 + c_4 k_4 +c_5 k_5 \\ y_{n+1}^* &= y_n + c_1^* k_1 + c_2^* k_2 +c_3^* k_3 + c_4^* k_4. \end{align} Because we can assume that the fifth-order solution is closer to the true solution than the fourth order algorithm, the difference $|y-y^*|$ provides a good estimate of the error of the fourth-order method. If this estimated error is larger than the desired tolerance, then the step size is decreased. If the error is smaller than the desired tolerance, the step size is increased. The Cash-Karp 5(4) ODE solver implements this technique for choosing the optimal step size.

In applications where the accuracy of the numerical solution is important, adaptive time step algorithms should always be used. As stated in *Numerical Recipes*:

*Many small steps should tiptoe through treacherous terrain, while a few great strides should speed through uninteresting countryside. The resulting gains in efficiency are not mere tens of percents or factors of two; they can sometimes be factors of ten, a hundred, or more.*

Adaptive step size algorithms are not well suited for tabulating functions or for simulation because the intervals between data points are not constant. An easy way to circumvent this problem is to use a method that takes multiple adaptive steps while checking to insure that the last step does not overshoot the desired fixed step size. ODE solvers often implement this technique. The solver acts like a fixed step size solver, even though the solver monitors its internal step size so as to achieve the desired accuracy.

It also is possible to combine the results from a calculation using two different values of the time step to yield a more accurate expression. Consider the Taylor series expansion of $f(t+dt)$ about $t$: $$\label{eq:motion/taylor1} f(t + dt) = f(t) + f'(t)dt + {1 \over 2!}f''(t) (dt)^2 + \ldots$$ Similarly, we have $$\label{eq:motion/taylor-1} f(t - dt) = f(t) - f'(t)dt + {1 \over 2!}f''(t) (dt)^2 + \ldots$$ We subtract (\ref{eq:motion/taylor-1}) from (\ref{eq:motion/taylor1}) to find the usual central difference approximation for the derivative $$\label{eq:motion/centraldifference} f'(t) \approx D_1(dt) = {f(t + dt) - f(t - dt) \over 2dt} - {(dt)^2 \over 6}f'''(t).$$ The truncation error is order $(dt)^2$. Next consider the same relation, but for a time step that is twice as large. $$\label{eq:motion/4more} f'(t) \approx D_1(2dt) = {f(t + 2dt) - f(t - 2dt) \over 4dt} - {4(dt)^2 \over 6}f'''(t).$$ Note that the truncation error is again order $(dt)^2$, but is four times bigger. We can eliminate this error to leading order by dividing (\ref{eq:motion/4more}) by 4 and subtracting it from (\ref{eq:motion/centraldifference}): \begin{align} f'(t) - {1 \over 4} f'(t) &= {3 \over 4} f'(t) \approx D_1(dt)- {1 \over 4} D_1(2dt), \nonumber \\ \mbox{or}\quad f'(t) &\approx {4D_1(dt) - D_1(2dt) \over 3}. \label{eq:motion/5point} \end{align} It is easy to show that the error for $f'(t)$ is order $(dt)^4$. Recursive difference formulas for derivatives can be obtained by canceling the truncation error at each order. This method is called *Richardson extrapolation*.

Another class of algorithms are *predictor-corrector* algorithms. The idea is to first *predict* the value of the new position: $$\label{eq:motion/xp} x_p = x_{n-1} + 2v_n dt. \qquad \qquad \mbox{(predictor)}$$ The predicted value of the position allows us to predict the acceleration $a_p$. Then using $a_p$, we obtain the *corrected* values of $v_{n+1}$ and $x_{n+1}$: \begin{align} v_{n+1}& = v_n + \frac{1}{2}(a_p + a_n) dt \\ x_{n+1}& = x_n + \frac{1}{2}(v_{n+1} + v_n) dt. \qquad \mbox{(corrected)} \end{align} The corrected values of $x_{n+1}$ and $v_{n+1}$ are used to obtain a new predicted value of $a_{n+1}$, and hence a new predicted value of $v_{n+1}$ and $x_{n+1}$. This process is repeated until the predicted and corrected values of $x_{n+1}$ differ by less than the desired value.

Note that the predictor-corrector algorithm is not self-starting. The predictor-corrector algorithm gives more accurate positions and velocities than the leapfrog algorithm and is suitable for very accurate calculations. However, it is computationally expensive, needs significant storage (the forces at the last two stages, and the coordinates and velocities at the last step), and becomes unstable for large time steps.

As we have emphasized, there is no single algorithm for solving Newton's equations of motion that is superior under all conditions. It is usually a good idea to start with a simple algorithm, and then to try a higher order algorithm to see if any real improvement is obtained.

We now discuss an important class of algorithms, known as *symplectic* algorithms, which are particularly suitable for doing long time simulations of Newton's equations of motion when the force is only a function of position. The basic idea of these algorithms derives from the Hamiltonian theory of classical mechanics. We first give some basic results needed from this theory to understand the importance of symplectic algorithms.

In Hamiltonian theory the generalized coordinates, $q_i$, and momenta, $p_i$, take the place of the usual positions and velocities familiar from Newtonian theory. The index $i$ labels both a particle and a component of the motion. For example, in a two particle system in two dimensions, $i$ would run from 1 to 4. The Hamiltonian (which for our purposes can be thought of as the total energy) is written as $$H({q_i,p_i}) = T + V, \label{eq:motion/hamiltonian}$$ where $T$ is the kinetic energy and $V$ is the potential energy. Hamilton's theory is most relevant for non-dissipative systems, which we consider here. For example, for a two particle system in two dimensions connected by a spring, $H$ would take the form: $$H = \frac{p_{1}^2}{2m} + \frac{p_{2}^2}{2m} + \frac{p_{3}^2}{2m} + \frac{p_{4}^2}{2m} + \frac{1}{2}k(q_{1} - q_{3})^2 + \frac{1}{2}k(q_{2} - q_{4})^2, \label{eq:motion/springH}$$ where if the particles are labeled as $A$ and $B$, we have $p_1 = p_{x,A}$, $p_2 = p_{y,A}$, $p_3 = p_{x,B}$, $p_4 = p_{y,B}$, and similarly for the $q_i$. The equations of motion are written as first-order differential equations known as Hamilton's equations: \label{eq:motion/hameq} \begin{align} {\dot p}_{i}& = -\frac{\partial H}{\partial q_i} \\ {\dot q}_{i}& = \frac{\partial H}{\partial p_i}, \end{align} which are equivalent to Newton's second law and an equation relating the velocity to the momentum. The beauty of Hamiltonian theory is that these equations are correct for other coordinate systems such as polar coordinates, and they also describe rotating systems where the momenta become angular momenta, and the position coordinates become angles.

Because the coordinates and momenta are treated on an equal footing, we can consider the properties of flow in phase space where the dimension of phase space includes both the coordinates and momenta. Thus, one particle moving in one dimension corresponds to a two-dimensional phase space. If we imagine a collection of initial conditions in phase space forming a volume in phase space, then one of the results of Hamiltonian theory is that this volume does not change as the system evolves. A slightly different result, called the *symplectic* property, is that the sum of the areas formed by the projection of the phase space volume onto the planes, ${q_i,p_i}$, for each pair of coordinates and momenta also does not change with time. Numerical algorithms that have this property are called symplectic algorithms. These algorithms are built from the following two statements which are repeated $M$ times for each time step. \begin{align} p_i^{(k+1)}& =p_i^{(k)} + a_{k} F_i^{(k)}\delta t \label{eq:motion/symplectic.a} \\ q_i^{(k+1)}& =q_i^{(k)} + b_{k} p_i^{(k+1)} \delta t, \label{eq:motion/symplectic.b} \end{align} where $F_i^{(k)} \equiv - \partial V(q_i^{(k)})/\partial q_i^{(k)}$. The label $k$ runs from 0 to $M-1$ and one time step is given by $\Delta t = M \delta t$. (We will see that $\delta t$ is the time step of an intermediate calculation that is made during the time step $dt$.) Note that in the update for $q_i$, the already updated $p_i$ is used. For simplicity, we assume that the mass is unity.

Different algorithms correspond to different values of $M$, $a_k$, and $b_k$. For example, $a_0 = b_0 = M = 1$ corresponds to the Euler-Cromer algorithm, and $M = 2$, $a_0 = a_1 = 1$, $b_0 = 2$, and $b_1 = 0$ is equivalent to the Verlet algorithm as we will now show. If we substitute in the appropriate values for $a_k$ and $b_k$ into Eqs. (\ref{eq:motion/symplectic.a}) and (\ref{eq:motion/symplectic.b}), we have \begin{align} p_i^{(1)}& = p_i^{(0)} + F_i^{(0)}\delta t \label{eq:motion/verletsym.a}\\ q_i^{(1)}& = q_i^{(0)} + 2p_i^{(1)}\delta t \label{eq:motion/verletsym.b}\\ p_i^{(2)}& = p_i^{(1)} + F_i^{(1)}\delta t \label{eq:motion/verletsym.c}\\ q_i^{(2)}& = q_i^{(1)}. \label{eq:motion/verletsym.d} \end{align} We next combine (\ref{eq:motion/verletsym.a}) and (\ref{eq:motion/verletsym.c}) for the momentum coordinate and (\ref{eq:motion/verletsym.b}) and (\ref{eq:motion/verletsym.d}) for the position, and obtain \begin{align} p_i^{(2)}&= p_i^{(0)} + (F_i^{(0)} + F_i^{(1)})\delta t \label{eq:motion/verletsym2.a}\\ q_i^{(2)}&= q_i^{(0)} + 2p_i^{(1)}\delta t. \label{eq:motion/verletsym2.b} \end{align} We take $\delta t = \Delta t/2$ and combine (\ref{eq:motion/verletsym2.b}) with (\ref{eq:motion/verletsym.a}) and find \begin{align} p_i^{(2)}&= p_i^{(0)} + \frac{1}{2}(F_i^{(0)} + F_i^{(1)})\Delta t \label{eq:motion/verletsym3.a}\\ q_i^{(2)}&= q_i^{(0)} +p_i^{(0)}\Delta t + \frac{1}{2} F_i^{(0)}(\Delta t)^2, \label{eq:motion/verletsym3.b} \end{align} which is identical to the Verlet algorithm, Eqs. (\ref{eq:motion/velverletx}) and (\ref{eq:motion/velverletv}), because for unit mass the force and acceleration are equal.

Reversing the order of the updates for the coordinates and the momenta also leads to symplectic algorithms: \begin{align} q_i^{(k+1)}& =q_i^{(k)} + b_{k}\delta t p_i^{(k)}, \label{eq:motion/symplectic2.a}\\ p_i^{(k+1)}& =p_i^{(k)} + a_{k}\delta t F_i^{(k+1)} \label{eq:motion/symplectic2.b} \end{align} A third variation uses different values of $k$ in one algorithm. Thus, if $M = 2$, which corresponds to two intermediate calculations per time step, we could use Eqs. (\ref{eq:motion/symplectic.a}) and (\ref{eq:motion/symplectic.b}) for the first intermediate calculation and Eqs. (\ref{eq:motion/symplectic2.a}) and (\ref{eq:motion/symplectic2.b}) for the second.

Why are these algorithms important? Because of the symplectic property, these algorithms will simulate an exact Hamiltonian, although not the one we started with in general. However, this Hamiltonian will be close to the one we wish to simulate if the $a_k$ and $b_k$ are properly chosen.

format: plain


text:
# Stats problems - Probability

Show 40 post(s) from this thread on one page
Page 2 of 2 First 12
• February 7th 2009, 02:05 AM
mr fantastic
Quote:

Originally Posted by mitch_nufc
Thanks, I done the corrections you suggested and i got exactly 1/3, would this be right?

Please show the exact calculation you did.

Here's the first line mine: $\frac{(0.95)(0.005)}{(0.95)(0.005) + (0.01)(0.995)}$.
• February 7th 2009, 02:11 AM
mitch_nufc
I eventually got what you got, my calculator has predictive fractions still on from when I was doing approximating fractions for irrational numbers, sorry about this
Show 40 post(s) from this thread on one page
Page 2 of 2 First 12


text:
17th International Conference on Topics in Astroparticle and Underground Physics (TAUP 2021)

26 August 2021 to 3 September 2021

Bridging the microhertz gap with asteroids: opportunities and challenges for gravitational wave detection

31 Aug 2021, 17:10
50m
Talk in parallel session Gravitational Waves

Speaker

Michael A. Fedderke (The Johns Hopkins University)

Description

The science case for a broad program of gravitational wave (GW) detection across all frequency bands is exceptionally strong. At present, there is a dearth of coverage by existing and proposed searches in the GW frequency band lying between the peak sensitivities of PTAs and LISA, roughly 0.1-100 microhertz. In this talk, I will outline a conceptual mission proposal to access this band. I will demonstrate that a few carefully chosen asteroids which orbit in the inner Solar System can act as excellent naturally occurring gravitational test masses despite the environmental noise sources. As such, a GW detector can be constructed by ranging between these asteroids using optical or radio links. At low frequencies, I will discuss how gravity gradient noise arising from the combined motion of the other $\sim 10^6$ asteroids in the inner Solar System sharply cuts off the sensitivity of this proposal. Sensitivity in the middle of this band is mostly limited by various solar perturbations to the asteroid test masses, while the high-frequency sensitivity is limited by noise in the ranging link. The projected strain-sensitivity curve that I will present indicates significant potential reach in this frequency band for a mission of this type.

Primary author

Michael A. Fedderke (The Johns Hopkins University)

Co-authors

Prof. Peter W. Graham (Stanford University) Prof. Surjeet Rajendran (The Johns Hopkins University)

Presentation Materials

 Fedderke_TAUP2021_GW_asteroids.pdf Recording of talk on Youtube
Your browser is out of date!

Update your browser to view this website correctly. Update my browser now

×


text:
# zbMATH — the first resource for mathematics

Defect correction methods for convection dominated convection-diffusion problems. (English) Zbl 0705.65081
The authors consider the approximate solution of singularly perturbed convection diffusion equations. Any appropriate numerical method faces the problem of resolution of sharp boundary layers. For instance, standard finite element methods are not appropriate when the perturbation parameter $$\epsilon$$ is of smaller order than the mesh distance h.
The method studied herein, originally due to P. W. Hemker [Lect. Notes Math. 960, 485-501 (1982; Zbl 0505.65047)], is a combination of defect correction with an artificial viscosity approximation. It computes a sequence of finite element solutions, where for $$\epsilon \ll h$$ the first solution is only a first order accurate approximation. At each correction step the residual is computed and a correction to the current approximation is calculated using the first order viscosity approximation.
The authors give local and global error estimates and show how the convergence depends on solution regularity, types of layers present and the subdomains on which the error is measured.
Reviewer: E.Lanckau

##### MSC:
 65N30 Finite element, Rayleigh-Ritz and Galerkin methods for boundary value problems involving PDEs 65N15 Error bounds for boundary value problems involving PDEs 35J70 Degenerate elliptic equations
Full Text:
##### References:
 [1] O. AXELSSON, On the numencal solution of convection dominated, convection-diffusion problems, in : Math. Meth. Energy Res. (K. I. Gross, ed. ), SIAM,Philadelphia, 1984. Zbl0551.76077 MR790509 · Zbl 0551.76077 [2] O. AXELSSON, Stability and error estimates of Galerkin finite element approximations for convection-diffusion equations, I. M. A. J. Numer. Anal., 1 (1981), 329-345. Zbl0508.76069 MR641313 · Zbl 0508.76069 · doi:10.1093/imanum/1.3.329 [3] W. ECKHAUS, Boundary layers in linear elliptic singular perturbation problems, SIAM Review, 14 (1972), 225-270. Zbl0234.35009 MR600325 · Zbl 0234.35009 · doi:10.1137/1014030 [4] V. ERVIN and W. LAYTON, High resolution minimal storage algorithms for convection dommated, convection diffusion equations, pp 1173-1201 in Tiams : of the Fourth Arms Conf. on Appl. Math. and Comp., 1987. Zbl0625.76095 MR905115 · Zbl 0625.76095 [5] V. ERVIN and W. LAYTON, An analysis of a defect correction method for a model convection diffusion equations, SIAM J. N. A. 26 (1989) 169-179. Zbl0672.65063 MR977954 · Zbl 0672.65063 · doi:10.1137/0726010 [6] P. W. HEMKER, Mixed defect correction iteration for the accurate solution of the convection diffusion equation, pp 485-501 in : Multigrid Methods, L. N. M. vol. 960, (W. Hackbusch and U. Trottenberg, eds.) Springer Verlag, Berlin 1982. Zbl0505.65047 MR685785 · Zbl 0505.65047 [7] P. W. HEMKER, The use of defect correction for the solution of a singularly perturbed o.d.e., preprint. CWI, Amsterdam, 1983. Zbl0504.65050 · Zbl 0504.65050 [8] C. JOHNSON and U. NÄVERT, An analysis of some finite element methods for advection diffusion problems, in : Anal. and Numer. Approaches to Asym. Probs. in Analysis (O. Axelson, L. S. Frank and A. van der Sluis, eds.) North Holland, 1981, 99-116. Zbl0455.76081 MR605502 · Zbl 0455.76081 [9] C. JOHNSON and U. NÄVERT and J. PITKARANTA, Finite element methods for linear hyperbolic problems, Comp. Meth. Appl. Mech. Eng., 45 (1984), 285-312. Zbl0526.76087 MR759811 · Zbl 0526.76087 · doi:10.1016/0045-7825(84)90158-0 [10] C. JOHNSON and A. H. SCHATZ and L. B. WAHLBIN, Crosswind smear and pointwise errors in streamline diffusion finite element methods, Math. Comp., 49 (1987), 25-38. Zbl0629.65111 MR890252 · Zbl 0629.65111 · doi:10.2307/2008248 [11] C. MIRANDA, Partial differential equations of elliptic type, Springer Verlag, Berlin, 1980. Zbl0198.14101 MR284700 · Zbl 0198.14101 [12] U. NÄVERT, A finite element method for convection diffusion problems, Ph. D. Thesis, Chalmers Inst. of Tech., 1982. [13] A. H. SCHATZ and L. WAHLBTN, On the finite element method for singularly perturbed reaction-diffusion problems in two and one dimensions, Math. Comp. 40 (1983), pp 47-89. Zbl0518.65080 MR679434 · Zbl 0518.65080 · doi:10.2307/2007363
This reference list is based on information provided by the publisher or from digital mathematics libraries. Its items are heuristically matched to zbMATH identifiers and may contain data conversion errors. It attempts to reflect the references listed in the original paper as accurately as possible without claiming the completeness or perfect precision of the matching.


text:
# LHC Shut Down By Weasel

It was reported last week that the LHC will be shut down for a few weeks due to a weasel gnawing through a power cable. The weasel, as you might expect, did not survive.

# Edelweiss Releases New WIMP Limits

The Edelweiss experiment has released a new preprint including updated limits in the low mass region from their WIMP dark matter search using germanium bolometers. There are still some stronger limits, but one interesting thing is that this further bolsters the case that the various purported WIMP signals from other experiments are probably not actually dark matter. The results from experiments like CoGeNT, CRESST, and DAMA all lie above the limit shown here.

# New Article on Sterile Neutrinos

Symmetry has a fairly new article on sterile neutrinos, explaining some of the basic ideas about what they are and why we’re interested in them. Measurements of the Z peak in e+/e- collisions at LEP showed that there are only 3 neutrinos, but there are some caveats to that. The measurements showing that there are 3 neutrinos really mean that there are 3 neutrinos that (1) have less than half the Z mass (91.2 GeV) and (2) interact with Standard Model particles via the weak force.

If we can instead add some “neutrinos” that don’t interact through the weak force, then there’s room for more neutrinos. One of the main ways people search for them is to find problems with the standard picture of 3-neutrino oscillations. If there are sterile neutrinos that mix with the three usual flavors (electron, muon, and tau), then maybe we can find evidence of neutrinos oscillating into sterile neutrinos. That is, regular neutrinos seemingly disappearing altogether rather than changing from one type to another.

# New e+/e- Collider Turns On

SuperKEKB, an upgraded version of the KEK B electron/positron collider at the KEK lab in Tsukuba has started up. There are no collisions yet, but they are starting to run beams around the main ring. The Belle II experiment will use he SuperKEKB facility to study B physics – that is, the physics of B mesons (and other particles involving b quarks). These are particularly interesting because the lifetimes of many of these particles decay weakly with long enough lifetimes to actually measure how far they travelled. B physics allows for high precision tests of various aspects of the Standard Model, including things like CP violation in the quark sector (i.e. matter/antimatter differences), particle spectroscopy (measuring the properties of the various kinds of composite particles), and searches for various kinds of new physics.

# ABCs of Particle Physics

Symmetry magazine has an ABCs of Particle Physics feature right now, which provides fun explanations of various topics in physics (with a focus on high energy physics). Everything is meant to be accessible to basically anyone. You can check it out here.

# How Does LIGO Point Back to a Source?

LIGO published the position of their gravitational wave on the sky when they announced that they had found something, so I was wondering how this could be done. While LIGO has more information that just timing, I decided to consider the case of reconstructing a position based just on the relative timing from two detectors.

I started with a coordinate system with the origin at the center of Earth (assumed to be a perfect sphere here). The z-axis points toward the north pole, the x-axis toward the equator at the Greenwich meridian and the y-axis 90 degrees east of the x-axis. I can then place the two detectors in this coordinate system (easiest using spherical coordinates). The time separation for a gravitational wave along some axis given by a unit vector u(θ,φ) is just

$\Delta t(\theta,\phi) = \frac{R_\oplus}{c}({\bf v_1}-{\bf v_2})\cdot {\bf u(\theta,\phi)}$

where v1 and v2 are the unit vectors pointing toward the two detectors, R is the radius of Earth, and c is the speed of light (the speed of the waves). If each detector can measure the time of the wave to some uncertainty σ and the relative timing uncertainty (typically from GPS) is negligible, I can define some likelihood that u is the correct vector as

$L(\theta,\phi) = e^{-\frac{1}{2}\chi^2(\theta,\phi)}$

where

$\chi^2(\theta,\phi) = \frac{(\Delta t_{\rm meas} - \Delta t(\theta,\phi))^2}{2\sigma^2}$

Using this function for some test values, we can try to understand how pointing works. If the time separation is 0, then we expect to see most likely values along the great circle between the two detectors.

I will also note that this is plotting position on the sky in Earth coordinates. Earth is not in an inertial reference frame: it rotates and moves in space, so actually mapping this against the stars in the sky gets very complicated and I’ll ignore that since it doesn’t affect what I want to show. The coordinates are also spherical coordinates, which are non-Euclidean, so shapes look strange if you’re not used to them. The curve above is actually the spherical version of a “straight” line. For a time separation of 8 ms, we see a thin ring that widens out as the resolution is degraded.

At 10 ms, we’re getting a large separation, and the ring has shrunk almost to just a blob.

Basically, what seems to be going on is:

1. The general shape of the likelihood distribution is a ring of positions corresponding to the measured time separation
2. As the time separation increases, the size of the ring on the sky decreases
3. When the separation is too large (requires a speed less than c), no position works well, so there is a sharp peak at the best (but still non-optimal) position.
4. The ring is characteristic of having two detectors. The center of the ring never changes (there are actually two centers depending on the sign of the separation) and its position is determined by the positions of the detectors
5. A third detector should break the degeneracy and allow for reconstruction of the position as a single point on the sky
6. Since the measured position looks like a blob, it looks like the time separation is large enough for the timing uncertainty to prevent us from seeing a ring-like shape.

# Finally Read the LIGO Paper

I finally found some time to read the LIGO paper. A couple things that I thought were interesting:

1. The peak power from gravitational waves was 200 solar masses per second. The power didn’t stay there for very long since a total of 3 solar masses was radiated away.
2. The rate of false positives the size of the signal seen is one in tens of thousands of years, so this is a signal that is enormously above any known backgrounds.
3. LIGO also uses some complicated template fitting routine where they compare the measured signal to a library of pre-calculated theoretical curves. This only gives approximate results for physics parameters, so they then have to supplement this with an actual fit.
4. The next biggest event had a false positive rate of only one every few years


text:
# Problem: Bismuth-210 is beta emitter with a half-life of 5.0 days. If a sample contains 1.3 g of Bi-210 (atomic mass = 209.984105 amu), how many beta emissions would occur in 13.6 days?

⚠Our tutors found the solution shown to be helpful for the problem you're searching for. We don't have the exact solution yet.

###### Problem Details

Bismuth-210 is beta emitter with a half-life of 5.0 days. If a sample contains 1.3 g of Bi-210 (atomic mass = 209.984105 amu), how many beta emissions would occur in 13.6 days?