Runge-Kutta 4 implementation in Rust

Question

This is a toy problem. I want to try using Rust for numerical modeling, but I found zero tutorials for this kind of stuff.

I also just started learning Rust yesterday, so please bear with me.

Usually when I'm doing an ODE solver, I create a grid array, then values array, define initial conditions and other parameters, then use a loop to find the solution.

Of course, I also need to write the result to a file, so I can use it later.

Now, I found it surprisingly hard to do any of these tasks in Rust. I eventually managed, but it took me most of the day to write this simple program.

Could you please check my code and point me to any possible improvements? Including the array creation, the loop, and file output.

// Program to find f(t) = cos(t) through ODE
// df/dt = g
// dg/dt = -f
use std::fs;
use std::io::Write;
fn main() {
 //make directory
 let dirname = "results";
 fs::create_dir(&dirname).expect("Error creating directory");
 //initiate the grid
 const NT:usize = 1000;
 let dt = 0.05;
 //initiate the time array
 let t: [f64; NT] = core::array::from_fn(|i| i as f64 * dt);
 //initiate exact solution array
 let fe: [f64; NT] = core::array::from_fn(|i| (i as f64 * dt).cos());
 //initiate numerical solution array
 //initial condition
 let f0 = 1.0;
 let g0 = 0.0;
 let mut f = [f0; NT];
 let mut g = [g0; NT];
 //for RK4 scheme
 let (mut k1, mut k2, mut k3, mut k4): (f64, f64, f64, f64);
 let (mut l1, mut l2, mut l3, mut l4): (f64, f64, f64, f64);
 for i in 1..NT {
 // RK4 scheme
 k1 = g[i-1];
 l1 = -f[i-1];
 k2 = k1 + l1*dt/2.0;
 l2 = l1 - k1*dt/2.0;
 k3 = k1 + l2*dt/2.0;
 l3 = l1 - k2*dt/2.0;
 k4 = k1 + l3*dt;
 l4 = l1 - k3*dt;
 // next values
 f[i] = -l1 + (k1 + 2.0*k2 + 2.0*k3 + k4)*dt/6.0;
 g[i] = k1 + (l1 + 2.0*l2 + 2.0*l3 + l4)*dt/6.0;
 }
 
 let filepath = format!("{}/{}",dirname,"RK4.dat");
 fs::File::create(&filepath).expect("Error creating file");
 let mut myfile = fs::OpenOptions::new()
 .write(true)
 .append(true)
 .open(filepath)
 .unwrap();
 let mut iline: String;
 for j in 0..NT {
 iline = format!("{} {} {}",t[j], fe[j], f[j]);
 writeln!(myfile, "{}", iline).expect("Error writing to file");
 }
}

This is the plot of RK4.dat, which shows that the program output is correct:

enter image description here

Answer 1

Use Rustfmt

Rust ships with a formatter, use it. It'll help keep your code readable and consistent with other code. Invoke it via cargo fmt or rustfmt path/to/file.rs.

Use Clippy

Rust also ships with Clippy, which contains a whole bunch lints to catch redundant, incorrect, or otherwise smelly code that doesn't make sense to have in rustc. I ran it on your code and got

warning: the borrowed expression implements the required traits
 --> src/main.rs:10:20
 |
10 | fs::create_dir(&dirname).expect("Error creating directory");
 | ^^^^^^^^ help: change this to: `dirname`
 |
 = help: for further information visit https://rust-lang.github.io/rust-clippy/master/index.html#needless_borrow
 = note: `#[warn(clippy::needless_borrow)]` on by default
warning: `playground` (bin "playground") generated 1 warning

Note that string literals are just an &'static str not a String, so they implement Copy, meaning that borrowing it is entirely useless here.

Avoid mixing core and std imports

For consistency's sake, unless you are writing code that needs to know about the distinction get all your imports from std, not core. In this case, It's using core::array::from_fn that you need to fix.

Use snake_case instead of runonsentencecase

Just a readability thing, rust uses snake case for variable names. For example, let myfile = /* stuff */ should be let my_file = /* stuff */.

File systems are hard

File systems are racy, complicated, platform-specific, and come with extremely awkward path handling. Rust does its best to provide tools to overcome these challenges, but they can be tricky to figure out for the uninitiated.

Use PathBuf to build paths

For any operation on paths more complicated than hardcoding, pull out std::path to handle the details for you. In particular, this will handle things like proper canonicalization and path separator differences for you. In your case that probably looks like this:

let mut filepath = PathBuf::from("results");
fs::create_dir(&filepath).expect("Error creating directory");
filepath.push("RK4.dat");
fs::File::create(&filepath).expect("Error creating file");
// and so on

Handle results directory already existing

fs::create_dir returns an error if the desired directory already exists, meaning your code unconditionally panics if run twice in the same directory. You can handle this by checking the returned error and continuing if it's caused by the directory already existing, or by using fs::create_dir_all, which considers the directory already existing being a success condition.

File::create doesn't do what you think it does

And even if it did, you've still opened yourself up to a race condition. File::create(path) works just like

let file = fs::OpenOptions::new()
 .write(true)
 .truncate(true)
 .create(true)
 .open(path);

This means your code creates/truncates and then opens the file, closes the file, and then immediately opens it in append mode without ensuring that it actually still exists. That isn't ideal- depending on whether you want to append or truncate the file either add .create(true) to your OpenOptions and drop File::create entirely or just use the value returned from File::create.

Inline your formatting when using writeln!

The writeln! macro supports the same formatting options that all rust's other formatting macros do, meaning the write loop should be:

for j in 0..NT {
 writeln!(myfile, "{} {} {}", t[j], fe[j], f[j]).expect("Error writing to file");
}

Flush writers before dropping them

Any time use you use an io::Write based object, you should flush it before dropping. This ensures data won't be silently lost if it can't be flushed successfully when it's being dropped.

// all your other file-based stuff here
myfile.flush().expect("failed to flush data to file");

Avoid panicking because of user error in application code

You should reserve panicking for instances when you the writer screw up, not for when you the user screw up. In your case, that means avoiding expect and unwrap calls near file system operations. I would recommend simply propagating the errors up through main, see the rust book's chapter on working with Results for a good explanation about that.

Consider buffering file writes

Not mandatory, but small numbers of large writes tend to be much faster when working with files then large numbers of small writes, so it may be wise to wrap the file in a BufWriter.

Final code and thoughts

// Program to find f(t) = cos(t) through ODE
// df/dt = g
// dg/dt = -f
use std::fs;
use std::io::{self, Write};
use std::path::PathBuf;
fn main() -> io::Result<()> {
 //initiate the grid
 const NT: usize = 1000;
 let dt = 0.05;
 //initiate the time array
 let t: [f64; NT] = std::array::from_fn(|i| i as f64 * dt);
 //initiate exact solution array
 let fe: [f64; NT] = std::array::from_fn(|i| (i as f64 * dt).cos());
 //initiate numerical solution array
 //initial condition
 let f0 = 1.0;
 let g0 = 0.0;
 let mut f = [f0; NT];
 let mut g = [g0; NT];
 //for RK4 scheme
 let (mut k1, mut k2, mut k3, mut k4): (f64, f64, f64, f64);
 let (mut l1, mut l2, mut l3, mut l4): (f64, f64, f64, f64);
 for i in 1..NT {
 // RK4 scheme
 k1 = g[i - 1];
 l1 = -f[i - 1];
 k2 = k1 + l1 * dt / 2.0;
 l2 = l1 - k1 * dt / 2.0;
 k3 = k1 + l2 * dt / 2.0;
 l3 = l1 - k2 * dt / 2.0;
 k4 = k1 + l3 * dt;
 l4 = l1 - k3 * dt;
 // next values
 f[i] = -l1 + (k1 + 2.0 * k2 + 2.0 * k3 + k4) * dt / 6.0;
 g[i] = k1 + (l1 + 2.0 * l2 + 2.0 * l3 + l4) * dt / 6.0;
 }
 let mut file_path = PathBuf::from("results");
 fs::create_dir_all(&file_path)?; //note that the `?` postfix operator is for error propagation, see the previously linked book chapter for details
 file_path.push("RK4.dat");
 
 let mut my_file = fs::File::create(&file_path)?;
 for j in 0..NT {
 writeln!(my_file, "{} {} {}", t[j], fe[j], f[j])?;
 }
 my_file.flush()?;
 Ok(())
}

Playground

Your code is pretty good overall, especially considering you've only been using rust for a day- the issues primarily come from the fact that file systems suck and code examples love to avoid proper error handling. Anyway, since you're interested in numerical modeling, here are a few crates that might be of interest:

• The num crate- a collection of utilities including big ints, complex and rational numbers, and abstractions over numerical operations and primitive types
• Plotters- a drawing library primarily designed for creating graphs, charts, and the like.
• Nalgebra- a general purpose linear algebra library. Also of note are rapier and salva, which are physics and fluid simulation libraries, respectively, and are created by the same organization.
• Ndarray- general purpose numpy style multidimensional array.
• Rayon- simple, safe, and fast parallel iterators. Not really specific to numerical modeling, but can greatly accelerate the processing of large quantities of data.

Anyway, welcome to rust, I hope you enjoy your stay.

• Yuriy S: Thank you! This is very helpful. I tried to fix all the issues you mentioned myself, before looking at your code. I kept getting an error trying to create a file through OpenOptions.new(), the only fix was to drop truncate(true), because apparently it conflicted with some other option ("Invalid argument" error). I see now that you just used File::create, does it mean I don't need to specify append() to be able to write lines to file?
• Aiden4: @YuriyS To write to a file it just has to be opened with write access. File::create does this automatically. Setting append(true) ensures the file is a) opened in write access and b} moves the file's cursor to the end so it doesn't overwrite old contents, whereas the cursor normally starts at the beginning of the file. This doesn't matter if you clear the file before use which is what File::create does and what I believe you want to do.