《四階龍格庫(kù)塔法解振蕩方程c源碼及畫圖》由會(huì)員分享,可在線閱讀,更多相關(guān)《四階龍格庫(kù)塔法解振蕩方程c源碼及畫圖(3頁(yè)珍藏版)》請(qǐng)?jiān)谘b配圖網(wǎng)上搜索。
1、
Vc++6.0 實(shí)現(xiàn)微分方程求解,程序如下:
#include
#include
double f1(double t, double x, double y, double z)
{
return -x + 6 * y + 2 * z;
}
double f2(double t, double x, double y, double z)
{
return -y + 3 * z + 2 * sin(t);
}
double f3(double t, double x, double y,
2、double z)
{
return -z + pow(t,2) * exp(-t) + cos(t);
}
double RungeKutta(double x0, double y0, double z0, double a, double b, int n, double x[], double y[], double z[], double (*f1)(double,double,double,double), double (*f2)(double,double,double,double), double (*f3)(double,double,double,do
3、uble))
{
double k1, k2, k3, k4 ;
double m1, m2, m3, m4 ;
double n1, n2, n3, n4 ;
double h;
double t;
int count;
count =0;
h = (b - a) / n;
x[0] = x0;
y[0] = y0;
z[0] = z0;
for(int i = 1;i<=n;i++)
{
t = a + (i - 1) * h;
4、 k1 = f1(t, x[i - 1], y[i - 1], z[i - 1]);
m1 = f2(t, x[i - 1], y[i - 1], z[i - 1]);
n1 = f3(t, x[i - 1], y[i - 1], z[i - 1]);
k2 = f1(t+h/2.0, x[i - 1] + k1 * h / 2.0, y[i - 1] + m1 * h / 2.0, z[i - 1] + n1 * h / 2.0);
m2 = f2(t+h/2.0, x[i - 1] + k1
5、* h / 2.0, y[i - 1] + m1 * h / 2.0, z[i - 1] + n1 * h / 2.0);
n2 = f3(t+h/2.0, x[i - 1] + k1 * h / 2.0, y[i - 1] + m1 * h / 2.0, z[i - 1] + n1 * h / 2.0);
k3 = f1(t+h/2.0, x[i - 1] + k2 * h / 2.0, y[i - 1] + m2 * h / 2.0, z[i - 1] + n2 * h / 2.0);
m3 = f2
6、(t+h/2.0, x[i - 1] + k2 * h / 2.0, y[i - 1] + m2 * h / 2.0, z[i - 1] + n2 * h / 2.0);
n3 = f3(t+h/2.0, x[i - 1] + k2 * h / 2.0, y[i - 1] + m2 * h / 2.0, z[i - 1] + n2 * h / 2.0);
k4 = f1(t+h/2.0, x[i - 1] + k3 * h / 2.0, y[i - 1] + m3 * h / 2.0, z[i - 1] + n3 * h / 2.
7、0);
m4 = f2(t+h/2.0, x[i - 1] + k3 * h / 2.0, y[i - 1] + m3 * h / 2.0, z[i - 1] + n3 * h / 2.0);
n4 = f3(t+h/2.0, x[i - 1] + k3 * h / 2.0, y[i - 1] + m3 * h / 2.0, z[i - 1] + n3 * h / 2.0);
x[i] = x[i - 1] + h * (k1 + 2 * k2 + 2 * k3 + k4) / 6.0;
y[i] =
8、y[i - 1] + h * (m1 + 2 * m2 + 2 * m3 + m4) / 6.0;
z[i] = z[i - 1] + h * (n1 + 2 * n2 + 2 * n3 + n4) / 6.0;
++count;
cout<<"t= "<