Не понимаю почему, мной написанный метод интегрирования Рунге-Кутта, считает не правильно. С каждым шагом погрешность увеличивается, до невозможного. Подскажите пожалуйста кто что знает...
Мой метод...
[codebox]#include "stdafx.h"
#include<iostream>
#include<math.h>
FILE *fout;
using namespace System;
using namespace std;
void print_y(double t, double y1, double y2, FILE* fout)
{
printf("t=%g\ty1=%g\ty2=%g\n",t,y1,y2);
fprintf(fout,"%f\t%f\t%f\n",t,y1,y2);
}
double f1(double y1, double y2)
{
return y2; //Ф-ия 1
}
double f2(double y1, double y2)
{
//return exp(y1)-exp(y2);
return -y1; //ф-ия 2
}
int main(array<System::String ^> ^args)
{
fout=fopen("data.txt","w");
fprintf(fout,"t\ty1\ty2\n");
double A, B; // пользовательский ввод
const int r=2;
double t;
double h;
int i;
cout<<"enter left bound: ";
cin>>A;
cout<<"enter right bound: ";
cin>>B;
Console::Write("Введите h= ");
cin>>h;
int m = (B-A)/h+1;
double** y = new double * [m];
for(i=0; i<m; i++)
y = new double[2];
cout<<"y1[0]= ";
cin>>y[0][0];
cout<<"y2[0]= ";
cin>>y[0][1];
t = A;
print_y(t,y[0][0],y[0][1],fout);
double k1,k2,k3,k4;
for(i=0;i<m-1;i++)
{
t+=h;
k1=f1(y[0]*h,y[1]);
k2=f1(y[0]+k1*h/4,y[1]+k1*h/4);
k3=f1(y[0]+h/10*(k1+4*k2),y[1]+h/10*(k1+4*k2));
y[i+1][0]=y[0]+h/3*(2*k1+4*k2+5*k3);
k1=f2(y[0]*h,y[1]);
k2=f2(y[0]+k1*h/4,y[1]+k1*h/4);
k3=f2(y[0]+h/10*(k1+4*k2),y[1]+h/10*(k1+4*k2));
y[i+1][1]=y[1]+h/3*(2*k1+4*k2+5*k3); //ур-е
print_y(t,y[i+1][0],y[i+1][1],fout);// печать y
}
return 0;
}[/codebox]
Мой метод...
[codebox]#include "stdafx.h"
#include<iostream>
#include<math.h>
FILE *fout;
using namespace System;
using namespace std;
void print_y(double t, double y1, double y2, FILE* fout)
{
printf("t=%g\ty1=%g\ty2=%g\n",t,y1,y2);
fprintf(fout,"%f\t%f\t%f\n",t,y1,y2);
}
double f1(double y1, double y2)
{
return y2; //Ф-ия 1
}
double f2(double y1, double y2)
{
//return exp(y1)-exp(y2);
return -y1; //ф-ия 2
}
int main(array<System::String ^> ^args)
{
fout=fopen("data.txt","w");
fprintf(fout,"t\ty1\ty2\n");
double A, B; // пользовательский ввод
const int r=2;
double t;
double h;
int i;
cout<<"enter left bound: ";
cin>>A;
cout<<"enter right bound: ";
cin>>B;
Console::Write("Введите h= ");
cin>>h;
int m = (B-A)/h+1;
double** y = new double * [m];
for(i=0; i<m; i++)
y = new double[2];
cout<<"y1[0]= ";
cin>>y[0][0];
cout<<"y2[0]= ";
cin>>y[0][1];
t = A;
print_y(t,y[0][0],y[0][1],fout);
double k1,k2,k3,k4;
for(i=0;i<m-1;i++)
{
t+=h;
k1=f1(y[0]*h,y[1]);
k2=f1(y[0]+k1*h/4,y[1]+k1*h/4);
k3=f1(y[0]+h/10*(k1+4*k2),y[1]+h/10*(k1+4*k2));
y[i+1][0]=y[0]+h/3*(2*k1+4*k2+5*k3);
k1=f2(y[0]*h,y[1]);
k2=f2(y[0]+k1*h/4,y[1]+k1*h/4);
k3=f2(y[0]+h/10*(k1+4*k2),y[1]+h/10*(k1+4*k2));
y[i+1][1]=y[1]+h/3*(2*k1+4*k2+5*k3); //ур-е
print_y(t,y[i+1][0],y[i+1][1],fout);// печать y
}
return 0;
}[/codebox]