#include <stdio.h>
#include <stdafx.h>
#include <math.h>
#include <iostream>
using namespace std;
double f(int i,double x,double y[4]){
switch (i){
case 1:return y[4];break;
case 2: return y[1];break;
case 3:return y[2];break;
case 4:return y[3];break;
default : break;
}
}
/*************************************************************************
Один шаг метода Рунге-Кутта четвертого порядка для решения
системы дифферециальных уравнений.
procedure SystemRungeKuttStep(
const x : Real;
const h : Real;
const n : Integer;
var y : array of Real);
Алгоритм совершает один шаг метода для системы
диффуров y'=F(i,x,y) для i=1..n
Начальная точка имеет кординаты (x,y[1], ..., y[n])
После выполнения алгоритма в переменной y содержится состояние
системы в точке x+h
*************************************************************************/
void step(double x,double h,int n,double y[4])
{
int i=0;
double yt[4];
double k1[4];
double k2[4];
double k3[4];
double k4[4];
for(i = 1; i <= n; i++)
{
k1 = h*f(i, x, y);
}
for(i = 1; i <= n; i++)
{
yt = y+0.5*k1;
}
for(i = 1; i <= n; i++)
{
k2 = h*f(i, x+h*0.5, yt);
}
for(i = 1; i <= n; i++)
{
yt = y+0.5*k2;
}
for(i = 1; i <= n; i++)
{
k3 = h*f(i, x+h*0.5, yt);
}
for(i = 1; i <= n; i++)
{
yt[i] = y[i]+k3[i];
}
for(i = 1; i <= n; i++)
{
k4[i] = h*f(i, x+h, yt);
}
for(i = 1; i <= n; i++)
{
y[i] = y[i]+(k1[i]+2.0*k2[i]+2.0*k3[i]+k4[i])/6;
}
}
/*************************************************************************
Алгоритм решает систему диффуров y[i]'=F(i,x,y) для i=1..n
методом Рунге-Кутта 4 порядка.
Начальная точка имеет кординаты (x,y[1], ..., y[n])
До конечной точки мы добираемся через n промежуточных
с постоянным шагом h=(x1-x)/m
Результат помещается в переменную result[4]
*************************************************************************/
void solvesystemrungekutta(double x,double x1,int steps,double result[4]){
for(int i = 1; i <= steps-1; i++)
{
step(x+i*(x1-x)/steps, (x1-x)/steps,4, result);
}
}
int _tmain(){
//первая координата
double temporaryresult[4];
solvesystemrungekutta(2,3,3,temporaryresult);
for(int i=1;i<=4;i++){
printf("Htpekmnf",temporaryresult[i],'\n');
}
char c=getchar();
return 0;
}
у меня вопрос как решить этим методом сисетму диффуров?
d2_x/dt_2=-GMx/(x^2+y^2)^(3/2)
d2_y/dy_2=-GMy/(x^2+y^2)^(3/2)