Runge-Kutta 4th order to solve 2nd order ODE

Feb 28, 2018 at 5:50am
I tried to write a code to solve a 2nd order ODE but somehow it did not work as I intended.

the equation is 2y" + 5y' +3y = 0 ,y(0) = 3 ,y'(0) = -4
The final answer will be y(x) = 3-4x
therefore y(1) = -1 and y'(1) = -4
But the program output is y(1) = -0.999995 which is somehow correct
However, y'(1) is output a wrong answer.
Please help!
Thank you!


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
#include <iostream>
#include <cmath>
#include <math.h> 

double ddx(double x,double y,double z)
{
	double dx = (-5*z-3*y)/2;
	return dx;
}

double test2ndorder(double x0, double y0, double z0, double x, double h)
{
	int n = (int)((x - x0) / h);
	double k1, k2, k3, k4;
	double l1, l2, l3, l4;

	double y = y0;
	double z = z0;

	for (int i = 1; i <= n; i++)
	{

		k1 = h * z0;
		l1 = h * ddx(x0, y0, z0);
		k2 = h * (z0 + 0.5*l1);
		l2 = h * ddx(x0 + 0.5*h, y0 + 0.5*k1, z0 + 0.5*l1);
		k3 = h * (z0 + 0.5*l2);
		l3 = h * ddx(x0 + 0.5*h, y0 + 0.5*k2, z0 + 0.5*l2);
		k4 = h * (z0 + l3);
		l4 = h * ddx(x0 + h, y0 + k3, z0 + l3);

		y = y + (1.0 / 6.0)*(k1 + 2 * k2 + 2 * k3 + k4);
		z = z + (1.0 / 6.0)*(l1 + 2 * l2 + 2 * l3 + l4);
		x = x + h;
	}
	return z;
}

int main()
{
	double x0, y0, z0, x, y, z,h;
	x0 = 0;
	x = 1;
	y0 = 3;
	z0 = -4;
	h =0.0001;

	z = test2ndorder(x0, y0, z0, x, h);
	std::cout << z;
}

Last edited on Mar 1, 2018 at 12:35am
Feb 28, 2018 at 6:57am
It won't compile if you start your function names with a 2.

ddx doesn't correspond to your given equation (you shouldn't divide by 2).

How do you know what dy/dx is at the end? - you don't output it. (EDIT - now you've edited your code and do output it ... but not y!)

Pass n, not h, or you will be subject to rounding error.

Please include headers.
Last edited on Feb 28, 2018 at 8:53am
Feb 28, 2018 at 7:25am
Sorry, the code I post before was an unupdate version.
Please look at the new version.
The output is 1.4614 but the answer should be -4 according to the equation.
Feb 28, 2018 at 8:07am
ddx doesn't correspond to the differential equation that you cite. Either your differential equation should have a 2 before the second derivative, or the result of this function shouldn't end by dividing by 2.

In neither case is the solution y=3-4x. It has a solution in decaying exponentials.

Your function should return both y and z(=dy/dx). You will need to return them through reference parameters.
Last edited on Feb 28, 2018 at 8:15am
Mar 1, 2018 at 6:00am
Thank you for the answer!

Turn out, I also use x0 ,y0 and z0 to calculate x, y and z in every step.
after I change that ,it works :D
Last edited on Mar 1, 2018 at 6:03am
Topic archived. No new replies allowed.