Object-oriented numerical calculation

Introduction

This article is the third day article of Meiji University Advent Calender 2018. This is the second post on the third day. It is a bicycle operation Advent Calender. Thank you.

Target

-** Implement numerical calculations in an object-oriented manner. ** ** -** Understand object orientation through numerical calculations. ** **

** (Note 1) The language used is Processing (almost Java). The difference from Java is that the drawing surface is substantial. ** ** ** (Note 2) There is no such thing as speeding up the calculation because it is implemented in an object-oriented manner. There is also a section that describes the advantages and disadvantages, so refer to that. ** **

Structure of this article

--Introduction of numerical calculation method (Euler method): --If you know the Euler method, skip it and ok --Object-oriented terminology: --If you know the class and interface, skip it and ok --Object-oriented numerical calculation: --Main part. Please read.

Introduction of numerical calculation method (Euler method)

Introducing the Euler method. The Euler method is one of the schemes (numerical calculation methods). First, I will introduce it in a general form, and then I will introduce what happens if it is a specific problem. Given the following initial value problem for ordinary differential equations.

\frac{dx}{dt} = f(x, t), \quad x(0 ) = x_0

If you can solve it by hand calculation in the form of $ f $, you can do it, but if you can not solve it, it is quick to rely on numerical calculation. The Euler method is equivalent to discretizing the above standing equation and creating and solving the following recurrence formula. Note that $ dt $ is the time difference width, which is a parameter of how finely the numerical calculation is performed.

x_{n + 1} = x_n + dt \cdot f(x_n, \ dt \cdot n)

You can find $ x_n $ sequentially by this recurrence formula.

A simple example is the Malthus model, an ordinary differential equation that represents population growth.

\frac{dx}{dt} = ax, \quad x(0) = x_0

This can be easily solved by hand calculation as $ x (t) = x_0 e ^ {at} $, but when performing numerical calculation by Euler method, the initial value is $ x_0 $ and the following recurrence Just think about the formula.

x_{n + 1} = x_n + dt \cdot (ax_n)

After that, if you solve this, you can get an answer that seems to be "somehow" close to the real solution.

Object-oriented term

Let's take a quick look at object orientation. Roughly speaking, object-oriented is the idea of seeing "things" and "concepts" as one "class". A "class" that requests implementation of necessary items is called an "interface".

I will explain using a common example. Suppose you want to implement an object-oriented "dog". What we want to ask a "dog" is, for example, barking. So the "dog interface" has a bark (barking) function.

However, even if you say that you bark in a bite, how to bark differs depending on the breed of dog. So, let's create a "poodle class" that inherits the "dog interface" and actually "implement" how to bark.

After that, I want to keep it as a pet, so I want to give each poodle a name. If that happens, you can give it a name. Individual poodles are called "dog class instances".

Implementing this idea in Java looks like this:

Dog.java


interface Dog {
    //Request implementation from inherited class
    void bark();
}

class Poodle implements Dog {
    String name;
    
    Poodle (String name) {
        this.name = name;
    }

    //Implement all interface functions
    void bark () {
        System.out.println("CanCan");
    }
}

//I'm creating an instance of Poodle as a Dog type.
Dog alto = new Poodle("Alto");
alto.bark(); //CanCan

Object-oriented numerical calculation

Now, the main subject is "object-oriented numerical calculation", but let's see how to think by taking the Euler method of the Malthus model as an example. First, the image looks like this. If you don't understand after reading this section, please return to the figure below. Screen Shot 2018-12-03 at 01.56.52.png

Overall picture

The entire numerical calculation is roughly divided into a "calculation part" and a "visualization part". The "calculation part" further includes "model setting" and "scheme setting". For example, in this case, the Malthus model is solved by the Euler method. We will look at the details below and explain using the ** bottom-up method **, which will eventually be summarized as a whole.

Malthus model in ODE

The Malthus model is ** one of the innumerable initial value problems of ordinary differential equations **. In other words, the Malthus model ** class ** is a ** implementation ** of the ** interface ** of the initial value problem of ordinary differential equations.

So, let's create an "ODE interface" as a big interface for the initial value problem of ordinary differential equations. Since ordinary differential equations are generally multiplied by $ \ frac {dx} {dt} = f (x, t), \ quad x (0) = x_0 $, each that implements an "ODE interface" Request the model to "implement $ f $".

Then, a concrete Malthus model is created by implementing the "ODE interface". Since the Malthus model is defined as $ f (x, t) = ax $ using the parameter $ a $ that represents the rate of population growth, the parameter $ a $ and the initial value are variables.

When implemented with the above in mind, it becomes as follows.

public interface ODE {
    public float getInitialVal ();
    public float f (float x, float t);
}


public class Malthus implements ODE {
    //Parameters
    private float a;
    private float initialVal;
  
    public Malthus (float a, float initialVal) {
        this.a = a;
        this.initialVal = initialVal;
    }
    
    //Implemented ODE interface request
    public float getInitialVal () {
        return this.initialVal;
    }

    public float f (float x, float t) {
        return this.a * x;
    }
}

Euler method in numerical scheme

The Euler method is ** one of many numerical schemes **. For example, other numerical schemes for ordinary differential equations include the Runge-Kutta method. This time, I chose the Euler method from among them. In object-oriented terms, the Euler method ** class ** is a ** implementation ** of the ** interface **, which is the numerical calculation scheme for ordinary differential equations.

Let's create a "Scheme interface" as an interface for the numerical calculation scheme of ordinary differential equations. Here, consider the functions that each scheme wants to have in common. Considering that the numerical calculation scheme of ordinary differential equations is to construct a recurrence formula, the following two requirements can be considered.

--calcNextX (ODE model, now $ x $, now $ t $) --A function that returns the value of $ x $ at the next time $ t + dt $ from the given model and the current $ x $ and $ t $. --getSolution (ODE model) --A function that finds an approximate solution at all discrete times by using calcNextX sequentially

After that, you can create an "Euler method class" that implements the "Scheme interface". With $ dt $ as the time difference, calcNextX can implement a recurrence formula represented by $ x_ {n + 1} = x_n + dt \ cdot f (x_n, n \ cdot dt) $.

Let's actually program the above.

public interface Scheme {
    public float calcNextX(ODE sys, float currentX, float currentT);
    public float[] getSolution (ODE sys);
}

public class Euler implements Scheme {
    private float dt;
    private int timeStep;
  
    public Euler (float dt, int timeStep) {
        this.dt = dt;
        this.timeStep = timeStep;
    }
    
    public float calcNextX(ODE sys, float currentX, float currentT) {
        return currentX + sys.f(currentX, currentT) * this.dt;
    }
  
    public float[] getSolution(ODE sys) {
        float[] xs = new float[this.timeStep];
        xs[0] = sys.getInitialVal();
    
        for (int i = 1; i < this.timeStep ; i++) {
             xs[i] = calcNextX(sys, xs[i - 1], (float)i * this.dt);
        }
    
        return xs;
    }
}

Integrate ODE and Scheme

Did you notice that you have created the "ODE interface" and the "Scheme interface" independently? In "Scheme", there was a function that takes "ODE" as an argument, but I don't care how "ODE" is implemented.

In order to handle "ODE" and "Scheme" collectively, it is necessary to create a large framework called "numerical calculation class" as another layer above. It is a class that has a role like an administrator.

The content is extremely simple, you just have to have one "ODE" and one "Scheme", and then have an approximate solution obtained using the given "ODE" and "Scheme".

The implementation is as follows.

public class NumericalAnalysis {
  private ODE sys;
  private Scheme scheme;
  private float[] solution;
  
  public NumericalAnalysis (ODE sys, Scheme scheme) {
    this.sys = sys;
    this.scheme = scheme;
    this.solution = this.getResult();
  }
  
  private float[] getResult() {
    return this.scheme.getSolution(this.sys);
  }
}

Visualization

Anyway, I want to visualize the numerical calculation result. Since it is a writing style peculiar to Processing from here, detailed explanation is omitted, but only the object-oriented part will be explained.

Since visualization is also a concept, create a "Graphic class". This is a class that interprets the results obtained by numerical calculation and draws.

It should be noted here that, from the standpoint of the "Graphic class", ** it only draws the given float array, not the numerical calculation result. ** In object-oriented programming, the numerical calculation part and the visualization part should be completely separated.

For the time being, how to put a pretty poor visualization class. All you have to do is draw an axis and draw the trajectory of the approximate solution. If you can use a graphic library, you should definitely use it.

public class Graphic {
    public void drawAxis () {
       background(255);
       stroke(0);
       // x-axis
       line(30, 30, 30, 470);
       // t-axis
       line(30, 470, 470, 470);
    }
  
    public void drawGraph (float[] solution) {
        drawAxis();
    
        float maxVal = max(solution);
        float minVal = min(solution);
        float dt = 440 / (float)solution.length;
    
        for (int i = 0; i < solution.length - 1; i++) {
            float now = map(solution[i], minVal, maxVal, 465, 35);
            float next = map(solution[i + 1], minVal, maxVal, 465, 35);
            line(dt * i + 30, now, dt * (i + 1) + 30, next);
        }
    }
}

Integration of numerical calculation and visualization

The "numerical calculation class" was a class that saved the numerical calculation results using the given "ODE" and "Scheme". Also, the "visualization class" was like drawing an array of floats. It is necessary to integrate these and create a class that passes the numerical calculation result obtained by the "numerical calculation class" to the "visualization class".

The "Overall class" plays this role. The contents are simple and connect the "numerical calculation class" and the "visualization class" so that the process from numerical calculation to visualization can be done by just calling one function (doAll function).

public class Overall {
    private NumericalAnalysis na;
    private Graphic g;
  
    public Overall (ODE sys, Scheme scheme) {
        this.na = new NumericalAnalysis (sys, scheme);
        this.g = new Graphic();
    }
  
    public void doAll () {
        g.drawGraph(this.na.getResult());
    }
}

Run

Now that it's complete, all you have to do is execute it. If you throw an instance of Malthus model and an instance of Euler method to Overall and execute it with doAll, the numerical calculation result will be visualized.

void setup () is a Processing-specific way of writing, but you can think of it as a main function.

void setup() {
    //Screen size
    size(500, 500);
    
    ODE sys = new Malthus(1.0, 1.0);
    Scheme scheme = new Euler(0.01, 300);
    Overall oa = new Overall (sys, scheme);
  
    oa.doAll();
}

Execution result Screen Shot 2018-12-03 at 04.19.26.png

It seems that the result like an exponential function is obtained properly. (Please forgive that there is no scale or label on the axis.)

merit and demerit

How was the object-oriented numerical calculation? That's right, claiming that it would be easier to write without object-orientation. It may look a little complicated because the program was constructed by superimposing abstraction on abstraction.

However, as a result of creating an interface that abstracts ODE and scheme in this way, it should be easier to create another model or scheme that is not the Malthus model or Euler method. For example, if you want to implement with the 4th-order Runge-Kutta method instead of the Euler method, you can already visualize by implementing the "Runge-Kutta class" from the "Scheme interface" and creating an instance of it ** It will end up. In other words, ** you can write another program just by replacing one part **.

Evolving topics

Other than ODE (PDE, etc.)

This time, we have focused on the initial value problem of ordinary differential equations and their numerical calculation schemes, but there are other partial differential equations and various numerical calculation methods that accompany them. However, generalization is not so difficult because a considerable framework can be explained by the same framework as this time.

Scheme to Singleton

In this article I wrote the scheme as an ordinary class. This means that you can create multiple instances of the same scheme in a program. That's bad. The scheme is just a method, so there shouldn't be more than one of the same.

There is an idea called Singleton as a way to solve it. Simply put, it's an idea to ensure that at most one instance of that class is created in your program. For more information, see [Design Patterns (link to Wikipedia)](https://ja.wikipedia.org/wiki/%E3%83%87%E3%82%B6%E3%82%A4%E3%83%B3% E3% 83% 91% E3% 82% BF% E3% 83% BC% E3% 83% B3_ (% E3% 82% BD% E3% 83% 95% E3% 83% 88% E3% 82% A6% E3 See% 82% A7% E3% 82% A2)).

It's an example that it should be a Singleton, and if you look at the linked design patterns, you'll know what design patterns are used in the program in this article and which ones should be used. ..

MVC (Model / View / Controller)

Who wondered why it was divided into a "numerical calculation class" and a "visualization class"? In the visualization section, I emphasized that the "visualization class" only translates and draws the float array, not the numerical calculation result.

Behind this word is the concept of MVC. MVC is an abbreviation for Model View Controller. This concept is easy to understand if you take a game as an example. Model is the game hardware and software, View is the game screen, and Controller is the controller. These are connected to each other through humans and signals. Based on the information displayed on the screen, it is input to the controller, the hardware interprets the input information, and the result is passed to the screen. They do not all share information at the same level as each other, but translate and exchange information. Although it is a delivery cycle, it is the outside such as humans and cables that connect them.

In this numerical calculation program, there was no Controller, but the numerical calculation class was Model and the visualization class was View. So he said that the visualization class only draws a float array.

Polymorphism

Polymorphism is translated as polymorphism in Japanese. This time, Scheme was prepared as an interface, and a framework was prepared in which the Euler method and Runge-Kutta method were implemented. In it, I created an Euler method instance as follows.

Scheme euler = new Euler(0.01, 1.0);

By doing this, I was able to create a variable that has the Euler method behavior even though it is a Scheme type. In this way, an idea that can be regarded as a concrete type in reality under an abstract type is called polymorphism. Polymorphism has the role of combining the cases of programs into one.

Encapsulation

Encapsulation is simply black boxing. From the outside, it is a concept that allows you to use the class without worrying about the complicated implementation method inside. For example, this time, using the Overall class that puts everything together, everything was executed with just one doAll function. The user can set various parameters of the Malthus model and Euler method and doAll to get the result without knowing the contents.

Summary

I created a program to perform numerical calculation of ordinary differential equations using object orientation. Using that program, I understood what object-oriented programming is.

Recommended Posts

Object-oriented numerical calculation
Conditional numerical calculation
Is WebAssembly really fast? [Numerical calculation]
Object-oriented summary
Object-oriented programming
Ruby / Rust linkage (3) Numerical calculation with FFI
[Java] Object-oriented
Ruby / Rust linkage (4) Numerical calculation with Rutie
Ruby / Rust linkage (5) Numerical calculation with Rutie ② Bezier