ode-solving/solvers.cpp at master · jamesfielder/ode-solving · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
#include "MFunction.hpp"
#include <fstream>
#include <iostream>
#include "MVector.hpp"
#include "solvers.hpp"

MVector ivp_solvers::euler(MVector &y, MFunction &f, int j){
	int i;
	double x, h;
	h = (b-a)/steps;
	x = a; // For printing, these are set to their initial values
	y = y; // and the loop started from 1 instead of 0 to get the file right
	if (j == 1){ // when printing
		ofstream output("output_euler.csv");
		output << x << "," << y[0] << "," << y[1] << endl;
		for(i = 1; i < steps+1; i++){
			x = a + i*h;
			y = y + h*f(x, y);
			output << x << "," << y[0] << "," << y[1] << endl;
		}
		output.close();
	}
	else { // not printing
		for(i = 1; i < steps+1; i++){
			x = a + i*h;
			y = y + h*f(x, y);
		}
	}
	// cout << y << endl; // optional print
	return y;
}
MVector ivp_solvers::midpoint(MVector &y, MFunction &f, int j){
	int i;
	double x, h;
	h = (b-a)/steps;
	x = a;
	y = y;
	if (j == 1){ //printing
		ofstream output("output_midpoint.csv");
		output << x << "," << y[0] << "," << y[1] << endl;
		for(i = 1; i < steps+1; i++){
			x = a + i*h;
			y = y + h*f(x + 0.5*h, y + 0.5*h*f(x,y));
			output << x << "," << y[0] << "," << y[1] << endl;
		}
		output.close();
	}
	else { // not printing
		for(i = 1; i < steps+1; i++){
			x = a + i*h;
			y = y + h*f(x + 0.5*h, y + 0.5*h*f(x,y));
		}
	}
	// cout << y << endl; // optional print
	return y;
}
MVector ivp_solvers::runge(MVector &y, MFunction &f, int j){
	int i;
	double x, h;
	MVector k1, k2, k3, k4;
	h = (b-a)/steps;
	x = a;
	y = y;
	k1 = h*f(x,y);
	k2 = h*f(x + 0.5*h, y + 0.5*k1);
	k3 = h*f(x + 0.5*h, y + 0.5*k2);
	k4 = h*f(x + h, y + k3);
	if (j == 1){
		ofstream output("output_runge.csv");
		output << x << "," << y[0] << "," << endl; //<< y[1] << endl;
		for(i = 1; i < steps+1; i++){
			x = a + i*h;
			k1 = h*f(x,y);
			k2 = h*f(x + 0.5*h, y + 0.5*k1);
			k3 = h*f(x + 0.5*h, y + 0.5*k2);
			k4 = h*f(x + h, y + k3);
			y = y + (1./6.)*(k1 + 2*k2 + 2*k3 + k4);
			output << x << "," << y[0] << endl;
		}
		output.close();
	}
	else {
		for(i = 1; i < steps+1; i++){
			x = a + i*h;
			k1 = h*f(x,y);
			k2 = h*f(x + 0.5*h, y + 0.5*k1);
			k3 = h*f(x + 0.5*h, y + 0.5*k2);
			k4 = h*f(x + h, y + k3);
			y = y + (1./6.)*(k1 + 2*k2 + 2*k3 + k4);
		}
	}
	// cout << y << endl; // option print statement
	return y;
}