C++：简单四阶 Runge-Kutta 计算器

Question

我正在制作一个计算器来求解我们在 class 中做的一些四阶 Runge-Kutta 方程式，虽然我能够让计算器工作并且运行，但它给我的价值观不太正确，我不明白为什么。

#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;
}

现在为了测试它我用了这个问题：

其中 h 是步长，因为它为每次迭代生成一个新的 Y 值，所以在第一次迭代后，一个错误会导致答案滚雪球。我认为它可能无法在小数点后计算足够的点数，但正如您所看到的，我使用了长双精度数但无济于事。使用 setprecision 命令不起作用，它只是在末尾添加许多零来从技术上满足该命令。

在这张照片中，我有一个 Runge-Kutta 计算器，我在左边找到了白色的，我自己的程序运行ning 在左边。正如您所看到的，第一次迭代是正确的，但之后他们变得越来越错误。

Answer 1

在该方法的每次迭代中，您都会更新 Y 近似值，但不会更新 X。尝试替换

K = Ks(PoN,Ans,A,B,X,Y,h);

和

K = Ks(PoN,Ans,A,B,i-h,Y,h);

（即将 X 更改为 i-h）。您可能还想为您的变量考虑更有意义的名称，并将您的函数分成更小的部分。这两者都可以使类似的东西更容易被发现。

C++：简单四阶 Runge-Kutta 计算器

C++: Simple Fourth Order Runge-Kutta calculator

c++

math

differential-equations

runge-kutta