1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63
|
#include <iostream>
#include <cmath>
double function(double t, double u[2][1])
{
double dudt[2][1];
dudt[0][0] = 3 * u[0][0] + 2 * u[1][0] - (2 * pow(t,2) + 1) * exp(2 * t);
dudt[1][0] = 4 * u[0][0] + u[1][0] + (pow(t,2) + 2 * t - 4) * exp(2 * t);
return dudt[2][1];
}
int main()
{
using namespace std;
double function(double t, double u[2][1]);
double h, t0, u1_0, u2_0;
int nos, i, j;
h = 0.1;
nos = 10;
t0 = 0.0;
u1_0 = 1.0;
u2_0 = 1.0;
double t[nos + 1];
double u[2][nos + 1];
t[0] = t0;
u[0][0] = u1_0;
u[1][0] = u2_0;
for(i = 0; i < nos; i++)
{
t[i + 1] = t[i] + h;
double u_n[2][1];
u_n[0][0] = u[0][i];
u_n[1][0] = u[1][i];
double k1[2], k2[2], k3[2], k4[2];
double u_npl[2][1];
for(j = 0; j < 2; j++)
{
k1[j] = h * function(t[i],u_n[j][0]);
k2[j] = h * function(t[i] + 0.5 * h, u_n[j][0] + 0.5 * k1[j]);
k3[j] = h * function(t[i] + 0.5 * h, u_n[j][0] + 0.5 * k2[j]);
k4[j] = h * function(t[i] + h, u_n[j][0] + k3[j]);
u_npl[j][0] = u_n[j][0] + ((k1[j] + 2 * k2[j] + 2 * k3[j] + k4[j]) / 6);
}
u[0][i + 1] = u_npl[0][0];
u[1][i + 1] = u_npl[1][0];
cout << "u[0][" << i + 1 << "] is: " << u[0][i + 1] << endl;
cout << "u[1][" << i + 1 << "] is: " << u[1][i + 1] << endl;
}
return 0;
}
|