C++

Forum

 
Why am I only returning 1?

ccwtree11 (44)
I've pasted my entire code (minus the header), but the focus is simply lines 130 onward in both files.. Basically I'm creating an Initial Value Problem solver, but I'm struggling to code the Explicit Runge-Kutta method. My problem file outputs the approximation to a text file, but I only get the value 1 for some reason. The commented code in the one_step_march for RK outputs the correct approximation, but this is a shortcut to the 4th order method, which isn't good because I need this to work on general order. So if anyone sees why I'm not getting incremental values, please let me know.
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
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
 
#include "ScalarIVPSolver.h"

ScalarIVPSolver::ScalarIVPSolver()
{
	G.resize(3);
}
void ScalarIVPSolver::set_function_G(int k, double(*fG)(vector<double> &))
{
	G[k] = fG;
}

ScalarIVPSolver::~ScalarIVPSolver()
{
	;
}

//−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−

TaylorSeriesMethod::TaylorSeriesMethod(int m)
{
	order = m;

	inversefactorial.resize(order + 1);
	inversefactorial[0] = 1;
	inversefactorial[1] = 1;
	double val = 1.0;
	for (int k = 2; k <= order; k++)
	{
		val *= ((double)k);
		inversefactorial[k] = 1.0 / val;
	}

	G.resize(order + 1);
}

double TaylorSeriesMethod::one_step_march(double h, double tau, double ytau)
{
	vector<double> x;
	x.push_back(tau);
	x.push_back(ytau);
	// Calculation of derivatives y^(k), k=1,...,order.
	for (int k = 1; k <= order; k++)
	{
		double tempval;
		tempval = (*G[k])(x);
		x.push_back(tempval);
	}

	// Calculation of Taylor’s polynomial coefficients.
	vector<double> polynomialcoeff(order + 1);
	for (int k = 0; k <= order; k++)
	{
		polynomialcoeff[k] = x[k + 1] * inversefactorial[k];
	}

	// Calculation of the polynomial. Use Horner’s Method.
	double ytauplush = 0.0; // Initialize.
	for (int k = order; k >= 0; k--)
	{
		ytauplush = ytauplush * h + polynomialcoeff[k];
	}
	return ytauplush;

}

//−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−

double EulerScheme::one_step_march(double h, double tau, double ytau)
{
	vector<double> x;
	x.push_back(tau); //why is this here?? just use G[0]??
	x.push_back(ytau);
	double ytauplush;

	// implement here
	ytauplush = ytau + h * (G[1])(x); //how does this work??

	return ytauplush;
}

//−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−

double BackwardEulerScheme::one_step_march(double h, double tau, double ytau)
{
	vector<double> x;
	x.push_back(tau);
	x.push_back(ytau);
	double ytauplush = ytau + h * (G[1])(x);
	double r = ytauplush; // initial guess
	x.push_back(r);
	double ytauplushP = h*G[2](x);

	// implement here
	for (int i = 1; i < 100000; i++)
	{
		r = r - ((ytauplush - r) / (ytauplushP - 1));
		if (abs(ytauplush - r) < 1 / 10000000000)
			break;
	}
	ytauplush = ytau + h * r;

	return ytauplush;
}

//−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−

double TrapezoidalScheme::one_step_march(double h, double tau, double ytau)
{
	vector<double> x;
	x.push_back(tau);
	x.push_back(ytau);
	double ytauplush = ytau + h * (G[1])(x);
	double r = ytauplush;
	x.push_back(r);
	double ytauplushP = (h / 2)*G[2](x);

	// implement here
	for (int i = 1; i < 100000; i++)
	{
		r = r - ((ytauplush - r) / (ytauplushP - 1));
		if (abs(ytauplush - r) < 1 / 10000000000)
			break;
	}
	ytauplush = ytau + (h / 2)*((G[1](x)) + r);

	return ytauplush;
}

//−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−

ExplicitRKScheme::ExplicitRKScheme(int s, const double *aa, const double *ww,
	const double *cc)
{
	stage = s;
	a = aa;
	w = ww;
	c = cc;
	k.resize(s);
}

double ExplicitRKScheme::one_step_march(double h, double tau, double ytau)
{
	double ytauplush, tauplush = tau, sum = 0, ak = 0, y = ytau;
	vector<double> x;
	x.resize(2);
	int n = 0;
	
	for (int i = 0; i < stage; i++)
	{
		tauplush = tau + c[i] * h;
		x[0] = tauplush;
		for (int j = 0; j < i; j++)
		{
			ak = ak + a[n] * k[j];
			n++;
		}
		y = ytau + h * ak;
		ak = 0;
		x[1] = y;
		k[i] = G[1](x);
		sum = sum + w[i] * k[i];
	}
	ytauplush = ytau + h * sum; //this is always 1...why???
	/*
	x.push_back(tau);
	x.push_back(ytau);
	k[0] = G[1](x);
	x[0] = tau + h / 2;
	x[1] = ytau + (h / 2)*k[0];
	k[1] = G[1](x);
	x[1] = ytau + (h / 2)*k[1];
	k[2] = G[1](x);
	x[0] = tau + h;
	x[1] = ytau + h*k[2];
	k[3] = G[1](x);
	ytauplush = ytau + (h / 6)*(k[0] + 2*k[1] + 2*k[2] + k[3]); //while this works
	*/
	return ytauplush;
}

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
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
 
#include "ScalarIVPSolver.h"

// Solve y’ = y, with y(0) = 1

// G0 is just the value of y(tau)
double G0(vector<double> &x)
{
	return x[1];
}

// This is nothing but the right hand side of the ODE
double G1(vector<double> &x)
{
	return x[1];
}

double G2(vector<double> &x)
{
	return x[2];
}

double G3(vector<double> &x)
{
	return x[3];
}

int main()
{
	int order = 3;

	ScalarIVPSolver *ts;
	ts = new TaylorSeriesMethod(order);

	// Set mathematical expression for calculation of y^(k), k=0,...,order
	ts->set_function_G(0, G0);
	ts->set_function_G(1, G1);
	ts->set_function_G(2, G2);
	ts->set_function_G(3, G3);

	double t0 = 0.0;
	double y0 = 1.0;
	double te = 2.0;
	int N = 20;
	double h = (te - t0) / N;
	vector<double> t(N + 1);
	vector<double> Y(N + 1);
	t[0] = t0;
	Y[0] = y0;
	for (int j = 1; j <= N; j++)
	{
		Y[j] = ts->one_step_march(h, t[j - 1], Y[j - 1]);
		t[j] = t[j - 1] + h;
	}

	delete ts;

	ofstream fileout;
	fileout.open("solts.txt");
	fileout << setprecision(10) << setiosflags(ios::scientific | ios::showpos);
	for (int j = 0; j <= N; j++)
	{
		fileout << t[j] << " " << Y[j] << " " << exp(t[j]) << endl;
	}
	fileout.close();

	// Run again with Euler Scheme (first order)
	ts = new EulerScheme;
	ts->set_function_G(0, G0);
	ts->set_function_G(1, G1);
	for (int j = 1; j <= N; j++)
	{
		Y[j] = ts->one_step_march(h, t[j - 1], Y[j - 1]);
		t[j] = t[j - 1] + h;
	}

	delete ts;

	fileout.open("soleuler.txt");
	fileout << setprecision(10) << setiosflags(ios::scientific | ios::showpos);
	for (int j = 0; j <= N; j++)
	{
		fileout << t[j] << " " << Y[j] << " " << exp(t[j]) << endl;
	}
	fileout.close();

	// Run again with Backward Euler Scheme (first order)
	ts = new BackwardEulerScheme;
	ts->set_function_G(0, G0);
	ts->set_function_G(1, G1);
	ts->set_function_G(2, G2);
	for (int j = 1; j <= N; j++)
	{
		Y[j] = ts->one_step_march(h, t[j - 1], Y[j - 1]);
		t[j] = t[j - 1] + h;
	}

	delete ts;

	fileout.open("solbackeuler.txt");
	fileout << setprecision(10) << setiosflags(ios::scientific | ios::showpos);
	for (int j = 0; j <= N; j++)
	{
		fileout << t[j] << " " << Y[j] << " " << exp(t[j]) << endl;
	}
	fileout.close();

	// Run again with Trapezoidal Scheme (second order)
	ts = new TrapezoidalScheme;
	ts->set_function_G(0, G0);
	ts->set_function_G(1, G1);
	ts->set_function_G(2, G2);
	for (int j = 1; j <= N; j++)
	{
		Y[j] = ts->one_step_march(h, t[j - 1], Y[j - 1]);
		t[j] = t[j - 1] + h;
	}

	delete ts;

	fileout.open("soltrapezoid.txt");
	fileout << setprecision(10) << setiosflags(ios::scientific | ios::showpos);
	for (int j = 0; j <= N; j++)
	{
		fileout << t[j] << " " << Y[j] << " " << exp(t[j]) << endl;
	}
	fileout.close();


	// Run again with Runge-Kutta (variable order)
	int stage = 4;
	double *weight = new double[stage] { (1/6),(1/3),(1/3),(1/6) }; // standard RK4 Butcher Table
	double *a = new double[6]{ 0.5,0,0.5,0,0,1 };
	double *c = new double[stage] {0, 0.5, 0.5, 1};
	ts = new ExplicitRKScheme(stage, a, weight, c);
	ts->set_function_G(0, G0);
	ts->set_function_G(1, G1);
	for (int j = 1; j <= N; j++)
	{
		Y[j] = ts->one_step_march(h, t[j - 1], Y[j - 1]);
		t[j] = t[j - 1] + h;
	}

	delete ts;
	delete c;
	delete a;
	delete weight;

	fileout.open("solrk.txt");
	fileout << setprecision(10) << setiosflags(ios::scientific | ios::showpos);
	for (int j = 0; j <= N; j++)
	{
		fileout << t[j] << " " << Y[j] << " " << exp(t[j]) << endl;
	}
	fileout.close();

	return 0;
}
Last edited on
Handy Andy (5051)
Hello ccwtree11,

I noticed that you started your question with a green check mark. This tells people that you have an answer to your question. Not sure if you can change the check mark, but it might be worth looking into.

One thing I did notice is that you open and close your output file several times. My first thought is that every time you open the file the file pointer is at the beginning of the file and overwrites anything that might be there. You might want to use: fileout.open("solbackeuler.txt", std::ios::app); to append the file at the end not overwrite the file.

Hope that helps,

Andy
Topic archived. No new replies allowed.