Solving 2nd order diff. equation using RK4 method

Problem related to 2nd order diff. equ. using RK4 method

Solve the following equation by using the Runge-Kutta fourth-order formula:

\[\frac{d^2y}{dx^2}+b\frac{dy}{dx}+cy=0\] \[y(0) = 1, \frac{dy}{dx}_{x=0}=0\] \[\delta =b^2 −4c\]

for

1. b > 0, $\delta$ < 0,
2. b < 0,$\delta$ < 0,
3. b > 0,$\delta$ = 0 and
4. b < 0, $\delta$ = 0.

Plot y(t) vs t for all the cases within a suitable range of t.

Answer

Code

#include <iostream>
#include <cmath>
#include <fstream>
using namespace std;
// We introduce z=dy/dt;
float f(float t, float y, float z, float b, float c) 
{
	return z;
}
// Defining the equation as in the equation with b and c as input
float g(float t, float y, float z, float b, float c) 
{
	return (-b * z - c * y);
}

float rk2(float t0, float y0, float z0, float tf, float h, float b, float c)
{
	float k1, k2, k3, k4, m1, m2, m3, m4, y, z, t;
	ofstream fileout;
	fileout.open("2ndorderdiff.dat");
	for (t = t0; t <= tf; t = t + h)
	{
		k1 = h * f(t, y0, z0, b, c);
		m1 = h * g(t, y0, z0, b, c);
		k2 = h * f(t + h / 2.0, y0 + k1 / 2.0, z0 + m1 / 2.0, b, c);
		m2 = h * g(t + h / 2.0, y0 + k1 / 2.0, z0 + m1 / 2.0, b, c);
		k3 = h * f(t + h / 2.0, y0 + k2 / 2.0, z0 + m2 / 2.0, b, c);
		m3 = h * g(t + h / 2.0, y0 + k2 / 2.0, z0 + m2 / 2.0, b, c);
		k4 = h * f(t + h, y0 + k3, z0 + m3, b, c);
		m4 = h * g(t + h, y0 + k3, z0 + m3, b, c);
		y = y0 + (k1 + 2 * k2 + 2 * k3 + k4) / 6.0;
		z = z0 + (m1 + 2 * m2 + 2 * m3 + m4) / 6.0;
		y0 = y;
		z0 = z;
		fileout << t << " " << y0 << "  " << z0 << endl;
	}
	return (0);
}

int main()
{
	float t0, y0, z0, tf, h, b, c, delta;
	cout << "Put the initial value of t: ";
	cin >> t0;
	cout << "Put the initial value of y at t=t0: ";
	cin >> y0;
	cout << "Put the initial value of z at t=t0: ";
	cin >> z0;
	cout << "Put the last value of t: ";
	cin >> tf;
	cout << "Put the step size: ";
	cin >> h;
	cout << "Input the values of b, delta: ";
	cin >> b >> delta;
	c = (b * b - delta) / 4.0;
	cout << "The value of c: " << c << endl;
	rk2(t0, y0, z0, tf, h, b, c);
	cout << "The calculations are done.";
	return 0;
}

Output

(1) b > 0, $\delta$ < 0

Put the initial value of t: 0
Put the initial value of y at t=t0: 1
Put the initial value of z at t=t0: 0
Put the last value of t: 10
Put the step size: 0.01
Input the values of b,c: 1
-1
The value of c: 0.5
The calculations are done.

(2) b < 0, $\delta$ < 0

Put the initial value of t: 0
Put the initial value of y at t=t0: 1
Put the initial value of z at t=t0: 0
Put the last value of t: 10
Put the step size: 0.01
Input the values of b,c: 1
-1
The value of c: 0.5
The calculations are done.

(3) b > 0, $\delta$ = 0

Put the initial value of t: 0
Put the initial value of y at t=t0: 1
Put the initial value of z at t=t0: 0
Put the last value of t: 10
Put the step size: 0.01
Input the values of b,c: 1
-1
The value of c: 0.5
The calculations are done.

(4) b < 0, $\delta$ = 0

Put the initial value of t: 0
Put the initial value of y at t=t0: 1
Put the initial value of z at t=t0: 0
Put the last value of t: 10
Put the step size: 0.01
Input the values of b,c: 1
-1
The value of c: 0.5
The calculations are done.