c++mathdifferential-equationsrunge-kutta
C++: Simple Fourth Order Runge-Kutta calculator
Im making a calculator to solve a few fourth order Runge-Kutta equations we were doing in class, and while I was able to get the calculator to work and run, the values it gives me are not quite right and I cannot figure out why.
#include <iostream>
#include <cmath>
#include <iomanip>
#include <cctype>
#include <string>
#include <windows.h>
#include <iostream>
#include <cstdlib>
#include <fstream>
#include <sstream>
#include <vector>
#include <numeric>
#include <algorithm>
using namespace std;
double PI = acos(-1.0);
long double Ks(int,int,long double,long double,long double,long double,long double);
int main ()
{
char answer;
do
int Ans,PoN;
long double X,Y,A,B,h,Fin,K;
cout << "Please enter X and Y used for problem [X Y]: ";
cin >> X >> Y;
cout << "What is the step size? ";
cin >> h;
cout << "What are we approximating to? Y(T) where T = ";
cin >> Fin;
cout << "Which equation form are you using (Please use number)?" << endl;
cout << "1) Ax (+/-) By" << endl;
cout << "2) Ax (+/-) By^2" << endl;
cout << "3) Ay (+/-) By^2" << endl;
cout << "Answer: ";
cin >> Ans;
cout << "You picked: ";
switch(Ans)
{
case 1:
cout << "Ax (+/-) By" << endl;
cout << "Enter A: ";
cin >> A;
cout << "Is the sign in the middle positive(1) or negative(2)? ";
cin >> PoN;
cout << "Enter B: ";
cin >> B;
break;
case 2:
cout << "Ax (+/-) By^2" << endl;
cout << "Enter A: ";
cin >> A;
cout << "Is the sign in the middle positive(1) or negative(2)? ";
cin >> PoN;
cout << "Enter B: ";
cin >> B;
break;
case 3:
cout << "Ay (+/-) By^2" << endl;
cout << "Enter A: ";
cin >> A;
cout << "Is the sign in the middle positive(1) or negative(2)? ";
cin >> PoN;
cout << "Enter B: ";
cin >> B;
break;
}
//cout << fixed << setprecision(8);
cout << " Y[#] | New Y | Exact" << endl;
cout << "------|---------|------" << endl;
for(long double i = X+h; i <= Fin; i += h)
{
K = Ks(PoN,Ans,A,B,X,Y,h);
Y = Y + (h/6)*K;
cout << setw(2) << "Y[" << setw(3) << i << setw(1) << "]" << "|" << Y << endl;
}
cout << "Would you like to repeat this program?(Y/N) ";
cin >> answer;
}
while (answer == 'y' || answer == 'Y');
if (answer == 'N' || answer == 'n')
{
cout << "Goodbye" << endl;
system("start notepad.exe");
}
system("pause");
return 0;
}
long double Ks(int PoN, int Ans, long double A,long double B,long double X,long double Y,long double h)
{
long double K1=0,K2=0,K3=0,K4=0,K=0;
switch(Ans)
{
case 1:
{
switch(PoN)
{
case 1: //Ax + By
K1 = A*X+B*Y;
K2 = A*(X+((1.0/2)*h)) + B*(Y+(1.0/2)*h*K1);
K3 = A*(X+((1.0/2)*h)) + B*(Y+(1.0/2)*h*K2);
K4 = A*(X+h) + B*(Y+h*K3);
break;
case 2: //Ax - By
K1 = (A*X)-(B*Y);
K2 = A*(X+((1.0/2)*h)) - B*(Y+(1.0/2)*h*K1);
K3 = A*(X+((1.0/2)*h)) - B*(Y+(1.0/2)*h*K2);
K4 = A*(X+h) - B*(Y+h*K3);
//cout << K1 << " " << K2 << " " << K3 << " " << K4 << endl;
break;
}
break;
}
case 2:
{
switch(PoN)
{
case 1: // Ax + By^2
K1 = A*X+pow(B*Y,2);
K2 = A*(X+((1.0/2)*h)) + pow(B*(Y+(1.0/2)*h*K1),2);
K3 = A*(X+((1.0/2)*h)) + pow(B*(Y+(1.0/2)*h*K2),2);
K4 = A*(X+h) + pow(B*(Y+h*K3),2);
break;
case 2: // Ax - By^2
K1 = A*X - pow(B*Y,2);
K2 = A*(X+((1.0/2)*h)) - pow(B*(Y+(1.0/2)*h*K1),2);
K3 = A*(X+((1.0/2)*h)) - pow(B*(Y+(1.0/2)*h*K2),2);
K4 = A*(X+h) - pow(B*(Y+h*K3),2);
break;
}
break;
}
case 3:
{
switch(PoN)
{
case 1: //Ay + By^2
K1 = A*X+B*Y;
K2 = A*(Y+((1.0/2)*h)) + B*(Y+(1.0/2)*h*K1);
K3 = A*(Y+((1.0/2)*h)) + B*(Y+(1.0/2)*h*K2);
K4 = A*(Y+h) + B*(Y+h*K3);
break;
case 2: //Ay - By^2
K1 = A*X-B*Y;
K2 = A*(Y+((1.0/2)*h)) - B*(Y+(1.0/2)*h*K1);
K3 = A*(Y+((1.0/2)*h)) - B*(Y+(1.0/2)*h*K2);
K4 = A*(Y+h) - B*(Y+h*K3);
break;
}
break;
}
}
K = (K1 + 2.0*K2 + 2.0*K3 + K4);
return K;
}
Now to test it I used this problem:
Where h is the step size, since it makes a new Y value for every iteration, one mistake causes the answer to snowball after the first iteration. I thought it might not be calculating enough points after the decimal but as you can see I used long doubles to no avail. Using the setprecision command does not work it just adds however many zeroes on the end to technically satisfy the command.
In this picture I have a Runge-Kutta calculator I found on the left in white and my own program running on the left. As you can see the first iteration is correct but after that they get more and more wrong.
Solution
With each iteration of the method, you update your Y approximation, but not the X. Try replacing
K = Ks(PoN,Ans,A,B,X,Y,h);
with
K = Ks(PoN,Ans,A,B,i-h,Y,h);
(i.e. change X to i-h).
You might also want to consider more meaningful names for your variables and dividing your functions into smaller pieces. Both of these could make something like this easier to spot.
- Fast & accurate atan/arctan approximation algorithm
- What's the difference between strtok_r and strtok_s in C?
- How memory address for pointer to arrays is same as an element in 2D array?
- Which is the best way to suppress "unused variable" warning
- How to use ellipsis in c's case statement?
- How can I exclude non-numeric keys? CS50 Caesar Pset2
- Fast ceiling of an integer division in C / C++
- Is there an invalid pthread_t id?
- How does SIMD (avx) processing work? for example, if I want 10 32 bit floats how do i fit in a 256 bit avx vector?
- FDCAN problems on STM32G4
- How does the call macro enable mutual recursion between functions f and g in this Hanoi Tower implementation?
- Running test on Rocket core CPU - global variable initialized to 0 is unsuccessful, output wrong value instead
- Interacting with C arrays without knowing the size
- Combination of two strings
- Avoiding strcpy overflow destination warning
- carriage return by fgets
- How to use special characters in C?
- Why does 1.0/100.0 == 0.1/10.0 give True?
- Is it correct to compare pointers in C?
- Force free() to return malloc memory back to OS
- How can I print to standard error in C with 'printf'?
- What is the standard behavior of fread in C on Windows?
- How is strtok removing lines it shouldn't have access to?
- Using array as smart point in C
- Assigning string to malloced 2d char array not working as intended
- How to refactor repetition inside a Makefile?
- Why does an empty preprocessor command still evaluate to something?
- How to implement variable sized array within C struct
- Character array typecasting to integer
- Handling HTTP Headers in a Minimal C HTTP Server