Есть задача.метод рунге.выдает значения.а как написать чтобы он еще график строил по этим значениям? помогите плиз))
Код | #include <conio.h> #include <math.h> #include <stdio.h>
typedef void (*ptr_func)(int, long double, long double*, long double*);
void rp1(int n, long double t, long double* y, long double* f) // Функция вычисления правой части { //коэффициенты long double u = 1; //средняя скорость течения потока (в день) long double k1 = 1; // long double k2 = 1; // long double k3 = 1; // long double db = 1; // long double la = 1; // long double cs = 1; // насыщенность (концентрация) растворенного кислорода
f[0] = 1.0/u *(-(k1 + k3) * y[0] + la); f[1] = 1.0/u *(-k1*y[0] + k2 * (cs - y[1]) - db); return; }
void rp2(int n, long double t, long double* y, long double* f) // Функция вычисления правой части { //коэффициенты long double k1 = 1; // long double k2 = 1; // long double k3 = 1; // long double db = 1; // long double la = 1; //
//f[0] = 1.0/u *(-(k1 + k3) * y[0] + la); //f[1] = 1.0/u *(-k1*y[0] + k2 * (cs - y[1]) - db);
//уравнения f[0] = (-k1 - k3) * y[0] + la; // биохимическая потребность растворенного кислорода в потоке f[1] = k1 * y[0] - k2 * y[1] + db; //концентрация растворенного кислорода на расстоянии вниз по течению потока
return; }
void method(ptr_func rp, int n, long double &t,long double &h, long double tk, long double eps, long double* y, long double r) // Метод Фельберга {
long double *yN, b1, bb1, bb2, t1, pp, ss; long double *k1, *k2, *k3, *k4, *k5, *k6; k1 = new long double [n]; // стадии метода k2 = new long double [n]; k3 = new long double [n]; k4 = new long double [n]; k5 = new long double [n]; k6 = new long double [n]; yN = new long double [n]; long double l1 =17./24.; long double l2 = 1./9.;
long double p51 = 16./135.; // коэффициенты метода long double p52 = 0.; long double p53 = 6656./12825.; long double p54 = 28561./56430.; long double p55 = -9./50.; long double p56 = 2./55.; long double p41 = 25./16.; long double p42 = 0.; long double p43 = 1408./2565.; long double p44 = 2197./4104.; long double p45 = -1./5.; long double p46 = 0;
long double b21 = 1./4.; long double b31 = 3./32.; long double b32 = 9./32.; long double b41 = 1932./2197.; long double b42 = -7200./2197.; long double b43 = 7296./2197.; long double b51 = 439./216.; long double b52 = -8.; long double b53 = 3680./513.; long double b54 = -845./4104.; long double b61 = -8./27.; long double b62 = 2.; long double b63 = -3544./2565.; long double b64 = 1859./4104.; long double b65 = -11./40.;
if(t+h > tk) h = fabs(tk - t);
rp(n, t, y, k1); do { t1 = t + b21*h; for(int i=0; i<n; i++) yN[i] = y[i] + b21*h*k1[i];
rp(n,t1,yN,k2);
t1 = t + (b31 + b32)*h; for (int i=0; i<n; i++) yN[i] = y[i] + h*(b31*k1[i] + b32*k2[i]); rp(n,t1,yN,k3);
t1 = t + (b41 + b42 + b43)*h; for (int i=0; i<n; i++) yN[i] = y[i] + h*(b41*k1[i] + b42*k2[i] + b43*k3[i]); rp(n,t1,yN,k4); t1 = t + (b51 + b52 + b53 + b54)*h; for(int i=0; i<n; i++) yN[i] = y[i] + h*(b51*k1[i] + b52*k2[i] + b53*k3[i] + b54*k4[i]); rp(n,t1,yN,k5);
t1 = t + (b61 + b62 + b63 + b64 + b65)*h; for(int i=0; i<n; i++) yN[i] = y[i] + h*(b61*k1[i] + b62*k2[i] + b63*k3[i] + b64*k4[i] + b65*k5[i]); rp(n,t1,yN,k6);
pp = 0.; ss = 0.; for(int i=0; i<n; i++) { b1 = fabs((p51 - p41)*k1[i] + (p52 - p42)*k2[i] + (p53 - p43)*k3[i] + (p54 - p44)*k4[i] + (p55 - p45)*k5[i] + (p56 - p46)*k6[i] )/(fabs(y[i]) + r); b1 = b1*l1; bb1 = fabs(k2[i] - k1[i]); bb2 = 0.e0; if(bb1 > 1.e-13) bb2 = fabs(32*k3[i] - 48*k2[i] + 16*k1[i])/bb1; if(pp < b1) pp = b1; if(ss < bb2) ss = bb2; } pp = h*pp; ss = ss*l2; if(pp < eps) break;
do { h = h/1.1; pp = pp/1.4641; if(pp <= eps) break;
h = h/1.1; } while (true); } while (true);
t=t + h; for(int i=0; i<n; i++) y[i] = y[i] + h*(p51*k1[i] + p52*k2[i] + p53*k3[i] + p54*k4[i] + p55*k5[i] + p56*k6[i]);
ss = 0.;
//m4: do { h = h*1.1;
pp = pp*1.61051; ss = ss*1.1; if (pp > eps || ss > 3.6) break; } while (true);
h = h/1.1;
delete yN; delete k1; delete k2; delete k3; delete k4; delete k5; delete k6; return; }
void main() { int n = 2; // размерность задачи
long double h=0.1; // начальный шаг интегрирования long double t= 0.0; // начальная точка интегрирования long double tk = 5; // конечная точка интегрирования long double eps = 0.01; // точность long double r = 1.; long double *y, *f; long double t1,h1;
y = new long double [n]; f = new long double [n];
y[0] = 1.e0; // начальные условия y[1] = 0.; printf("\tt\t h\t\ty1(t)\t y2(t)\n");
do { t1=t; method(rp1, n,t,h,tk,eps,y,r); h1 = t - t1; printf("\n%12.5f %12.5f %12.5f %12.5f",t,h1,y[0],y[1]); } while (fabs(t-tk)>1.e-13);
printf("\n\n");
//система (2) h=0.1; // начальный шаг интегрирования t= 0.0; // начальная точка интегрирования tk = 5; // конечная точка интегрирования eps = 0.01; // точность r = 1.; y[0] = 1.e0; // начальные условия y[1] = 0.; printf("\tt\t h\t\ty1(t)\t y2(t)\n");
do { t1=t; method(rp2, n,t,h,tk,eps,y,r); h1 = t - t1; printf("\n%12.5f %12.5f %12.5f %12.5f",t,h1,y[0],y[1]); } while (fabs(t-tk)>1.e-13);
getch(); } |
Модератор: не забываем пользоваться кнопочкой "Код". |