C ++: Простой калькулятор Рунге-Кутта четвертого порядка - PullRequest
0 голосов
/ 30 ноября 2018

Я создаю калькулятор для решения нескольких уравнений Рунге-Кутты четвертого порядка, которые мы делали в классе, и хотя я смог заставить калькулятор работать и работать, значения, которые он мне дает, не совсем верны, и я не могу понятьпочему

#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 не работает, она просто добавляет столько нулей на конце, чтобы технически удовлетворить команду.enter image description here

На этом снимке у меня калькулятор Рунге-Кутты, который я нашел слева белым, и моя собственная программа, работающая слева.Как видите, первая итерация верна, но после этого они все больше и больше ошибаются.

1 Ответ

0 голосов
/ 01 декабря 2018

С каждой итерацией метода вы обновляете свое Y приближение, но не X.Попробуйте заменить

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

на

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

(т. Е. Изменить X на i-h).Возможно, вы также захотите рассмотреть более значимые имена для ваших переменных и разделить ваши функции на более мелкие части.И то, и другое может облегчить обнаружение чего-то подобного.

...