* Данная работа не является научным трудом, не является выпускной квалификационной работой и представляет собой результат обработки, структурирования и форматирования собранной информации, предназначенной для использования в качестве источника материала при самостоятельной подготовки учебных работ.
“Реализация примера решений дифференциальног о уравнения второго пор ядка методом Р унга-Кутта при использовании компилятора C ++”
ВВЕДЕНИЕ
Численное интегрирование обыкновенных уравнений
В данной работе рассматривается уравне ние для задач с частым изменением шага , с учетом , что испол ьзуемый метод Ру нге-Кутта требует относительно большого количеств а вычислений производных на каждом шаге и для них весьма сложен эффективный контро ль величины шага.
При решении обыкновенных дифференциальных уравнений высших порядков необходимо учитыв ать , что каждое уравнение равносильно системе уравнений первого порядка . В частност и в данном случае при использовании метод а Рунге-Кутта рассматривается система обыкновенны х дифференциальных уравнений первого порядка вида
y ’ = f ( x , y , z , … ..)
с решением :
y = y ( x )
…
z = z ( x )
C учетом , что любая разностная схема из уравнений втор ого порядка может быть применена к каждом у из уравнений первого порядк а с обозначением :
y ( x k )= y k и тп.
f ( x k , y k , z k , … ..)= f k и тп .
Ввиду практической важности дифференциальных уравнений второго порядка п редставляют интерес схемы численного инте грирования дифференциального уравнения
вида :
y ” = f ( x , y , y ’ )
при наличии начальных условий :
x (0)
y (0)
y ’ (0)= V y (0)
x ’ (0)= V x (0)
В работе испльзовался только численный метод для уравнений второго порядка , без применения схем “ предсказание-коррекция” и интерполяционно-итерационной .
y k +1 = y k + y ’ k ∙'95∆ t+1/6(k 1 +k 2 +k 3 ) ∙'95∆ t y ’ k +1 = y ’ k + 1/6(k 1 +2k 2 +2k 3 + k 4 )
x k +1 = x k + x ’ k ∙'95∆ t +1/6( k 1 + k 2 + k 3 ) ∙'95∆ t x ’ k +1 = x ’ k + 1/6(k 1 +2k 2 +2k 3 + k 4 )
где :
k 1 = f ( x k , y k )
k 2 = f ( x k + V xk ∙'95∆ t /2, y k + V yk ∙'95∆ t /2)
k 3 = f ( x k + V xk ∙'95∆ t /2+( k 1 ∙'95∆ t )/2, y k + V yk ∙'95∆ t /2+( k 1 ∙'95∆ t )/2)
k 4 = f ( x k + V xk ∙'95∆ t +( k 2 ∙'95∆ t )/2, y k + V yk ∙'95∆ t +( k 2 ∙'95 ∆ t )/2)
ПРАКТИЧЕСКОЕ ПРИМЕНЕНИЕ
В работ е исп ользовался gcc , g ++ - GNU project C and C ++ Compiler ( v 2.7), FreeBSD radius . local . stv . ee 3.5- STABLE F reeBSD 3.5- STABLE , @ radius . local . stv . ee :/ usr / src / sys / compile / cmg i 386 .
В расчетах указывались постоянные : Y (гравитационная постоянная )=6,67*10 -11 Нм 2 /кг 2 ,
Dt =100 сек ., М (Земля )=5,796*10 24 кг при формуле :
YM /( x 2 + y 2 ) 3/2 y
YM /( x 2 + y 2 ) 3/2 x
Так при вводе вс е х значений =0, видим , что координаты не меняю тся , из чего следует , что тело при знач ениях =0 находится в состоянии покоя :
Enter value for x:0
Enter value for y:0
Enter value for Vx:0
Enter value for Vy:0
| N | X | Y | V x | Vy |
| 0 | x= 0 | y= 0 | Vx= 0 | Vy= 0 |
| 1 | x= 0 | y= 0 | Vx= 0 | Vy= 0 |
| 2 | x= 0 | y= 0 | Vx= 0 | Vy= 0 |
| 3 | x= 0 | y= 0 | Vx= 0 | Vy= 0 |
| 4 | x= 0 | y= 0 | Vx= 0 | Vy= 0 |
| 5 | x= 0 | y= 0 | Vx= 0 | Vy= 0 |
| 6 | x= 0 | y= 0 | Vx= 0 | Vy= 0 |
| 7 | x= 0 | y= 0 | Vx= 0 | Vy= 0 |
| 8 | x= 0 | y= 0 | Vx= 0 | Vy= 0 |
| 9 | x = 0 | y = 0 | Vx = 0 | Vy = 0 |
Далее приведем изменение значений при заданных параметрах , отличных от нуля :
Enter value for x:1
Enter value for y:2
Enter value for Vx:3
Enter value for Vy:4
| N | X | Y | Vx | Vy |
| 0 | x= 1 | y= 2 | Vx= 3 | Vy= 4 |
| 1 | x= -2147483648 | y= 402 | Vx= -1073 741821 | Vy= 4 |
| 2 | x= -2147483648 | y= 802 | Vx= -2147483648 | Vy= 4 |
| 3 | x= -2147483648 | y= -2147483648 | Vx= -2147483648 | Vy= 177585892 |
| 4 | x= -2147483648 | y= -1568763632 | Vx= -214 7483648 | Vy= 177585892 |
| 5 | x= -2147483648 | y= -2147483648 | Vx= -2147483648 | Vy= -811648420 |
| 6 | x= -2147483648 | y= -1707947024 | Vx= -2147483648 | Vy= -811648420 |
| 7 | x= -2147483648 | y= -2147483648 | Vx= -2147483648 | Vy= -1287506 838 |
| 8 | x= -2147483648 | y= -2049148568 | Vx= -2147483648 | Vy= -1287506838 |
| 9 | x= -2147483648 | y= -2147483648 | Vx= -2147483648 | Vy= 255514688 |
| 10 | x= -2147483648 | y= 1929148672 | Vx= -2147483648 | Vy= 255514688 |
PROGRAMM LISTING
// Koduto"o"u"lesanne
// author- Marina Mitrofanova
// kood- 951464
#include
#include
#include
#include
#define STEPS 20
#define M_M 5.97e24*6.67e-11
#define dt 100
#define STEPS 20
long calc_f(long c_p,long v_c_p,double k1_p,double k2_p,double k3_p)
return c_p+v_c_p*dt+(k1_p+k2_p+k3_p)*dt/6;
long prefix(long y_p,long x_p)
if(x_p==0||y_p==0)
return 0;
else
return M_M/pow((x_p*x_p+y_p*y_p),3/2);
int main()
long x,y,vx,vy;
double yk1=0,yk2=0,yk3=0,yk4=0;
double yk1_r=0,yk2_r=0,yk3_r=0,yk4_r=0;
double xk1=0,xk2=0,xk3=0,xk4=0;
double xk1_r=0,xk2_r=0,xk3_r=0,xk4_r=0;
int i;
printf("Enter value for x:");
scanf("%d", &x);
printf("Enter value for y:");
scanf("%d", &y);
printf("Enter v alue for Vx:");
scanf("%d", &vx);
printf("Enter value for Vy:");
scanf("%d", &vy);
printf("| N | X | Y | Vx | Vy |
for(i=0;i<=STEPS;i++)
printf("| %2i | x=%12d | y=%12d | Vx=%12d | Vy=%12d |\n",i,
yk1_r=prefix(x,y)*y;
yk2_r=prefix(x,y)*calc_f((y+vy*dt/2),vy,yk1,yk2,yk3);
yk3_r=prefix(x,y)*calc_f((y+vy*dt/2+yk1_r*dt/4),vy,yk1,yk2,y
yk4_r=prefix(x,y)*calc_f((y+vy*dt+yk2_r*dt/2),vy,yk1,yk2,yk3
yk1=yk1_r;yk2=yk2_r;yk3=yk3_r;yk4=yk4_r;
y=calc_f(y,vy ,yk1,yk2,yk3);
vy=vy+(yk1+2*yk2+2*yk3+yk4)/6;
xk1_r=prefix(x,y)*x;
xk2_r=prefix(x,y)*calc_f((x+vx*dt/2),vx,xk1,xk2,xk3);
xk3_r=prefix(x,y)*calc_f((x+vx*dt/2+xk1_r*dt/4),vx,xk1,xk2,x
xk4_r=prefix(x,y)*calc_f((x+vx*dt+xk2_r*dt/2),vx,xk1,xk2,xk3
xk1=xk1_r;xk2=xk2_r;xk3=xk3_r;xk4=xk4_r;
x=calc_f(x,vx,xk1,xk2,xk3);
vx=vx+(xk1+2*xk2+2*xk3+xk4)/6;