[C++] runge kutty felberg — błąd kompilacji
: 27 sie 2016, o 22:02
Wydzielenie zmiennych jest czytelniejsze dla osób mniej biegłych. Tu nie przekazywałem dobrych praktyk, a jedynie bezbłędny (z punktu widzenia kompilatora) kod.
Forum matematyczne: miliony postów, setki tysięcy tematów, dziesiątki tysięcy użytkowników - pomożemy rozwiązać każde zadanie z matematyki
https://matematyka.pl/
makosia15 pisze:Tak dobrze . Równanie różniczkowe drugiego rzędu rozbija się na równanie pierwszego rzędu. Stosuje się w takim przypadku układ równań. Różniczka drugiego rzędu, to jeden układ równań itd. Do scałkowania równań różniczkowych zwyczajnych stosuje się między innymi metodę rungego kutty felberga, który jest bardziej rozbudowany od metody runge kutty, ale działa na podobnej zasadzie. Są to metody całkowania różniczki zupełnej. W tych metodach wystarczy wstawić do równania podstawowego układ równań jako funkcje, algorytm też ulega w takim przypadku drobnym zmianom. W tym ostatnim linku co podałam to jest metoda runge kutty dla równania różniczkowego drugiego rzędu. Więc rozbita jest na układ równań, algorytm podstawowy też ulega drobnym zmianom. Nie wiem dokładnie, ale są tam tablice \(\displaystyle{ k}\) i \(\displaystyle{ k[j]}\), które odpowiadają pierwszej pochodnej \(\displaystyle{ k[1][j]}\) i drugiej pochodnej \(\displaystyle{ k[2][j]}\).
Kod: Zaznacz cały
int main()
Kod: Zaznacz cały
double f(double x, double y2)
{
return y2;
}
double f1(double x, double y1, double y2)
{
return (x+y_1+y_2)*(x+y_1+y_2);
}
Kod: Zaznacz cały
#include<iostream>
#include<cmath>
#include<fstream>
using namespace std;
double f1(double x, double y0, double y1)
{
return y1;
}
double f(double x, double y0, double y1)
{
return (x+y1+y0)*(x+y1+y0);
}
int main()
{
int n=10;
double x[n],y1[n],y0[n],y3[n],y2[n],k[3][7],c0[n],c1[n];
double h=0.1;
int i;
y0[0]=0;
y1[0]=0;
x[0]=0;
for (i=1;i<=n;i++)
{
x[i]=x[0]+i*h;
k[1][1]=h*f1(x[i],y0[i],y1[i]);
k[2][1]=h*f(x[i],y0[i],y1[i]);
k[1][2] = h* f1(x[i] + (h/4), y0[i] + (k[1][1]/4), y1[i]+(k[2][1]/4));
k[2][2] = h* f(x[i] + (h/4), y0[i] + (k[1][1]/4), y1[i]+(k[2][1]/4));
k[1][3]= h * f1(x[i] + (3*h/8), y0[i] + (3*k[1][1]/32) + (9*k[1][2]/32),y1[i] + (3*k[2][1]/32) + (9*k[2][2]/32));
k[2][3]= h * f(x[i] + (3*h/8), y0[i] + (3*k[1][1]/32) + (9*k[1][2]/32),y1[i] + (3*k[2][1]/32) + (9*k[2][2]/32));
k[1][4]= h * f1(x[i] + (12*h/13), y0[i] + (1932*k[1][1]/2197) - (7200*k[1][2]/2197) + (7296*k[1][3]/2197),y1[i] + (1932*k[2][1]/2197) - (7200*k[2][2]/2197) + (7296*k[2][3]/2197));
k[2][4]= h * f(x[i] + (12*h/13), y0[i] + (1932*k[1][1]/2197) - (7200*k[1][2]/2197) + (7296*k[1][3]/2197),y1[i] + (1932*k[2][1]/2197) - (7200*k[2][2]/2197) + (7296*k[2][3]/2197));
k[1][5]= h * f1(x[i] + h, y0[i] + (439*k[1][1]/216) - 8*k[1][2] + (3680*k[1][3]/513) - (845*k[1][4]/4104),y1[i] + (439*k[2][1]/216) - 8*k[2][2] + (3680*k[2][3]/513) - (845*k[2][4]/4104));
k[2][5]= h * f(x[i] + h, y0[i] + (439*k[1][1]/216) - 8*k[1][2] + (3680*k[1][3]/513) - (845*k[1][4]/4104),y1[i] + (439*k[2][1]/216) - 8*k[2][2] + (3680*k[2][3]/513) - (845*k[2][4]/4104));
k[1][6]= h * f1(x[i] + (h/2), y0[i] - (8*k[1][1]/27) + 2*k[1][2] - (3544*k[1][3]/2565) + (1859*k[1][4]/4104) - (11*k[1][5]/40),y1[i] - (8*k[2][1]/27) + 2*k[2][2] - (3544*k[2][3]/2565) + (1859*k[2][4]/4104) - (11*k[2][5]/40));
k[2][6]= h * f(x[i] + (h/2), y0[i] - (8*k[1][1]/27) + 2*k[1][2] - (3544*k[1][3]/2565) + (1859*k[1][4]/4104) - (11*k[1][5]/40),y1[i] - (8*k[2][1]/27) + 2*k[2][2] - (3544*k[2][3]/2565) + (1859*k[2][4]/4104) - (11*k[2][5]/40));
//cout<<k[1][1]<<" "<<k[1][2]<<" "<<k[1][3]<<" "<<k[1][4]<<" "<<k[1][5]<<" "<<k[1][6]<<endl;
//cout<<k[2][1]<<" "<<k[2][2]<<" "<<k[2][3]<<" "<<k[2][4]<<" "<<k[2][5]<<" "<<k[2][6]<<endl;
y3[i+1]=y1[i]+((25*(k[1][1]/216))+(1408*(k[1][3]/2565))+(2197*(k[1][4]/4101))-(k[1][5]/5));
y2[i+1]=y0[i]+((25*(k[2][1]/216))+(1408*(k[2][3]/2565))+(2197*(k[2][4]/4101))-(k[2][5]/5));
y1[i+1]=y1[i]+((16*(k[1][1]/135))+(6656*(k[1][3]/12852))+(28561*(k[1][4]/56430))-(9*(k[1][5]/50))+(2*(k[1][6]/55)));
y0[i+1]=y0[i]+((16*(k[2][1]/135))+(6656*(k[2][3]/12852))+(28561*(k[2][4]/56430))-(9*(k[2][5]/50))+(2*(k[2][6]/55)));
//c1[i+1] = y1[i+1] - y3[i+1];
cout << endl << x[i] << " " <<y1[i]<<" "<<y0[i]<< " "<<c1[i=1]<<endl << endl;
}
system("pause");
return 0;
}
Kod: Zaznacz cały
#include<iostream>
#include<cmath>
#include<fstream>
using namespace std;
double f1(double x, double y0, double y1)
{
return y1;
}
double f(double x, double y0, double y1)
{
return (x+y1+y0)*(x+y1+y0);
}
int main()
{
int n=10;
double x[n],y1[n],y0[n],y3[n],y2[n],k[3][7],c0[n],c1[n];
double h=0.1;
int i;
y0[0]=0;
y1[0]=0;
x[0]=0;
for (i=1;i<=n;i++)
{
x[i]=x[0]+i*h;
k[1][1]=h*f1(x[i],y0[i],y1[i]);
k[2][1]=h*f(x[i],y0[i],y1[i]);
k[1][2] = h* f1(x[i] + (h/4), y0[i] + (k[1][1]/4), y1[i]+(k[2][1]/4));
k[2][2] = h* f(x[i] + (h/4), y0[i] + (k[1][1]/4), y1[i]+(k[2][1]/4));
k[1][3]= h * f1(x[i] + (3*h/8), y0[i] + (3*k[1][1]/32) + (9*k[1][2]/32),y1[i] + (3*k[2][1]/32) + (9*k[2][2]/32));
k[2][3]= h * f(x[i] + (3*h/8), y0[i] + (3*k[1][1]/32) + (9*k[1][2]/32),y1[i] + (3*k[2][1]/32) + (9*k[2][2]/32));
k[1][4]= h * f1(x[i] + (12*h/13), y0[i] + (1932*k[1][1]/2197) - (7200*k[1][2]/2197) + (7296*k[1][3]/2197),y1[i] + (1932*k[2][1]/2197) - (7200*k[2][2]/2197) + (7296*k[2][3]/2197));
k[2][4]= h * f(x[i] + (12*h/13), y0[i] + (1932*k[1][1]/2197) - (7200*k[1][2]/2197) + (7296*k[1][3]/2197),y1[i] + (1932*k[2][1]/2197) - (7200*k[2][2]/2197) + (7296*k[2][3]/2197));
k[1][5]= h * f1(x[i] + h, y0[i] + (439*k[1][1]/216) - 8*k[1][2] + (3680*k[1][3]/513) - (845*k[1][4]/4104),y1[i] + (439*k[2][1]/216) - 8*k[2][2] + (3680*k[2][3]/513) - (845*k[2][4]/4104));
k[2][5]= h * f(x[i] + h, y0[i] + (439*k[1][1]/216) - 8*k[1][2] + (3680*k[1][3]/513) - (845*k[1][4]/4104),y1[i] + (439*k[2][1]/216) - 8*k[2][2] + (3680*k[2][3]/513) - (845*k[2][4]/4104));
k[1][6]= h * f1(x[i] + (h/2), y0[i] - (8*k[1][1]/27) + 2*k[1][2] - (3544*k[1][3]/2565) + (1859*k[1][4]/4104) - (11*k[1][5]/40),y1[i] - (8*k[2][1]/27) + 2*k[2][2] - (3544*k[2][3]/2565) + (1859*k[2][4]/4104) - (11*k[2][5]/40));
k[2][6]= h * f(x[i] + (h/2), y0[i] - (8*k[1][1]/27) + 2*k[1][2] - (3544*k[1][3]/2565) + (1859*k[1][4]/4104) - (11*k[1][5]/40),y1[i] - (8*k[2][1]/27) + 2*k[2][2] - (3544*k[2][3]/2565) + (1859*k[2][4]/4104) - (11*k[2][5]/40));
//cout<<k[1][1]<<" "<<k[1][2]<<" "<<k[1][3]<<" "<<k[1][4]<<" "<<k[1][5]<<" "<<k[1][6]<<endl;
//cout<<k[2][1]<<" "<<k[2][2]<<" "<<k[2][3]<<" "<<k[2][4]<<" "<<k[2][5]<<" "<<k[2][6]<<endl;
y3[i+1]=y1[i]+((25*(k[1][1]/216))+(1408*(k[1][3]/2565))+(2197*(k[1][4]/4101))-(k[1][5]/5));
y2[i+1]=y0[i]+((25*(k[2][1]/216))+(1408*(k[2][3]/2565))+(2197*(k[2][4]/4101))-(k[2][5]/5));
y1[i+1]=y1[i]+((16*(k[1][1]/135))+(6656*(k[1][3]/12852))+(28561*(k[1][4]/56430))-(9*(k[1][5]/50))+(2*(k[1][6]/55)));
y0[i+1]=y0[i]+((16*(k[2][1]/135))+(6656*(k[2][3]/12852))+(28561*(k[2][4]/56430))-(9*(k[2][5]/50))+(2*(k[2][6]/55)));
//c1[i+1] = y1[i+1] - y3[i+1];
cout << endl << x[i] << " " <<y1[i]<<" "<<y0[i]<< " "<<c1[i=1]<<endl << endl;
}
system("pause");
return 0;
}