用runge kutta方法求解颂歌系统 [英] Using runge kutta method to solve a system of odes

问题描述

大家好。



我的工作是在Runge Kutta方法和射击方法的帮助下解决常微分方程组。



以下是程序和结果窗口。当我按下调试时，结果只显示从7.0到10.0的x，但x的整个范围是从0.0到10.0：

似乎它没有计算任何东西。非常感谢你的帮助。



ODE系统：



f'''+ 3 * f * f'' - 2 * f'* f'+ q = 0

q''+ 2.25 * f * q'= 0

原始和边界条件：f （0）= f'（0）= 1-q（0）= f'（无穷大）= q（无穷大）= 0



结果窗口：

x [i] = 7.050000然后y [i] = - nan< ind>

x [i] = 7.050000然后z [i] = - nan< ind>

x [i] = 7.050000然后w [i] = - nan< ind>

x [i] = 7.050000然后p [i] = - nan< ind>

..................................等

< br $> b $ b

Hello everybody

My work is to solve a system of Ordinary Differential Equations with the help of Runge Kutta method and Shooting method.

Below is the program and result window. when I press debugging, the result only displays x from 7.0 to 10.0, however the full range of x is from 0.0 to 10.0:
It seems like It doesn't calculate anything. I really appreciate your help.

System of ODEs:

f'''+ 3*f*f''-2*f'*f'+q=0
q''+2.25*f*q'=0
Innitial and boundary conditions: f(0)=f'(0)=1-q(0)=f'(infinity)=q(infinity)=0

The result window:
"x[i]=7.050000 then y[i]=-nan<ind>
x[i]=7.050000 then z[i]=-nan<ind>
x[i]=7.050000 then w[i]=-nan<ind>
x[i]=7.050000 then p[i]=-nan<ind>
..................................etc"

#include<math.h>
double f1(double e1, double e2, double e3, double e4, double e5, double e6)
{
return e3;
}
//---------------------
double f2(double e1, double e2, double e3, double e4, double e5, double e6)
{
return e4;
}
//---------------------
double f3(double e1, double e2, double e3, double e4, double e5, double e6)
{
return -3 * e2*e4 + 2 * e3*e3 - e5;
}
//---------------------
double f4(double e1, double e2, double e3, double e4, double e5, double e6)
{
return e6;
}
double f5(double e1, double e2, double e3, double e4, double e5, double e6)
{
return -2.25*e2*e6;
}
#include<stdio.h>
#include<conio.h>
#include<iostream>
using namespace std;

//------------------------
int main(void)
{
const double EPS =1e-5;
int i, n = 200;
double h;
printf("number of interval in the coordinate n= %d\n", n);
//
//
//
h = 10.0 / n;
//
double *x;
x = new double[n + 1];
x[0] = 0.0;
//
for (i = 0; i<n + 1; i++)
x[i] = i*h;
//
//
printf("step in the coordinate h= %f\n", h);
//
double *u;
u = new double[n + 1];
double *y;
y = new double[n + 1];
y[0] = 0.0;
double *z;
z = new double[n + 1];
z[0] = 7.0;
double *w;
w = new double[n + 1];
w[0] = 1.0;
double *p;
p = new double[n + 1];
p[0] = -3.0;
//
int iter;
//
double eta = 7.0, epsilon = -3.0;
double hz, hp, dd, dz, dp;
double bb, cc, bbz, ccz, bbp, ccp;
double a1, a2, a3, a4;
double b1, b2, b3, b4;
double c1, c2, c3, c4;
double d1, d2, d3, d4;
double g1, g2, g3, g4;
iter = 0;
do
{
iter++;
u[0] = 0.0;
y[0] = 0.0;
z[0] = eta;
w[0] = 1.0;
p[0] = epsilon;
//
//
for (i = 0; i<n; i++)
{
a1 = f1(x[i], u[i], y[i], z[i], w[i], p[i])*h;
b1 = f2(x[i], u[i], y[i], z[i], w[i], p[i])*h;
c1 = f3(x[i], u[i], y[i], z[i], w[i], p[i])*h;
d1 = f4(x[i], u[i], y[i], z[i], w[i], p[i])*h;
g1 = f5(x[i], u[i], y[i], z[i], w[i], p[i])*h;
//
a2 = f1(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1)*h;
b2 = f2(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1)*h;
c2 = f3(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1)*h;
d2 = f4(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1)*h;
g2 = f5(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1)*h;
//
//
a3 = f1(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2)*h;
b3 = f2(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2)*h;
c3 = f3(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2)*h;
d3 = f4(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2)*h;
g3 = f5(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2)*h;
//
//
a4 = f1(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3)*h;
b4 = f2(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3)*h;
c4 = f3(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3)*h;
d4 = f4(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3)*h;
g4 = f5(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3)*h;

//
u[i + 1] = u[i] + (a1 + 2.0*a2 + 2.0*a3 + a4) / 6.0;
y[i + 1] = y[i] + (b1 + 2.0*b2 + 2.0*b3 + b4) / 6.0;
z[i + 1] = z[i] + (c1 + 2.0*c2 + 2.0*c3 + c4) / 6.0;
w[i + 1] = w[i] + (d1 + 2.0*d2 + 2.0*d3 + d4) / 6.0;
p[i + 1] = p[i] + (g1 + 2.0*g2 + 2.0*g3 + g4) / 6.0;
}
//
bb = y[n] - 0.0;
cc = w[n] - 0.0;
cout << "bb= " << bb << endl;
cout << "cc= " << cc << endl;
//
//
u[0] = 0.0;
y[0] = 0.0;
if (fabs(eta) <= 1.0) hz = sqrt(EPS);
if (fabs(eta)>1.0) hz = eta*sqrt(EPS);
z[0] = eta + hz;
w[0] = 1.0;
p[0] = epsilon;
//
//
for (i = 0; i<n; i++)
{
a1 = f1(x[i], u[i], y[i], z[i], w[i], p[i])*h;
b1 = f2(x[i], u[i], y[i], z[i], w[i], p[i])*h;
c1 = f3(x[i], u[i], y[i], z[i], w[i], p[i])*h;
d1 = f4(x[i], u[i], y[i], z[i], w[i], p[i])*h;
g1 = f5(x[i], u[i], y[i], z[i], w[i], p[i])*h;
//
//
a2 = f1(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1)*h;
b2 = f2(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1)*h;
c2 = f3(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1)*h;
d2 = f4(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1)*h;
g2 = f5(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1)*h;
//
a3 = f1(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2)*h;
b3 = f2(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2)*h;
c3 = f3(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2)*h;
d3 = f4(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2)*h;
g3 = f5(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2)*h;
//
a4 = f1(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3)*h;
b4 = f2(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3)*h;
c4 = f3(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3)*h;
d4 = f4(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3)*h;
g4 = f5(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3)*h;
//
u[i + 1] = u[i] + (a1 + 2.0*a2 + 2.0*a3 + a4) / 6.0;
y[i + 1] = y[i] + (b1 + 2.0*b2 + 2.0*b3 + b4) / 6.0;
z[i + 1] = z[i] + (c1 + 2.0*c2 + 2.0*c3 + c4) / 6.0;
w[i + 1] = w[i] + (d1 + 2.0*d2 + 2.0*d3 + d4) / 6.0;
p[i + 1] = p[i] + (g1 + 2.0*g2 + 2.0*g3 + g4) / 6.0;

}
bbz = y[n] - 0.0;
ccz = w[n] - 0.0;
//
//
u[0] = 0.0;
y[0] = 0.0;
z[0] = eta;
w[0] = 1.0;
if (fabs(epsilon) <= 1.0) hp = sqrt(EPS);
if (fabs(epsilon)>1.0) hp = epsilon*sqrt(EPS);
p[0] = epsilon + hp;
//
//
for (i = 0; i<n; i++)
{
a1 = f1(x[i], u[i], y[i], z[i], w[i], p[i])*h;
b1 = f2(x[i], u[i], y[i], z[i], w[i], p[i])*h;
c1 = f3(x[i], u[i], y[i], z[i], w[i], p[i])*h;
d1 = f4(x[i], u[i], y[i], z[i], w[i], p[i])*h;
g1 = f5(x[i], u[i], y[i], z[i], w[i], p[i])*h;
//
a2 = f1(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1)*h;
b2 = f2(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1)*h;
c2 = f3(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1)*h;
d2 = f4(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1)*h;
g2 = f5(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1)*h;
//
a3 = f1(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2)*h;
b3 = f2(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2)*h;
c3 = f3(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2)*h;
d3 = f4(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2)*h;
g3 = f5(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2)*h;
//
a4 = f1(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3)*h;
b4 = f2(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3)*h;
c4 = f3(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3)*h;
d4 = f4(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3)*h;
g4 = f5(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3)*h;
//
u[i + 1] = u[i] + (a1 + 2.0*a2 + 2.0*a3 + a4) / 6.0;
y[i + 1] = y[i] + (b1 + 2.0*b2 + 2.0*b3 + b4) / 6.0;
z[i + 1] = z[i] + (c1 + 2.0*c2 + 2.0*c3 + c4) / 6.0;
w[i + 1] = w[i] + (d1 + 2.0*d2 + 2.0*d3 + d4) / 6.0;
p[i + 1] = p[i] + (g1 + 2.0*g2 + 2.0*g3 + g4) / 6.0;


}
//
//
bbp = y[n] - 0.0;
ccp = w[n] - 0.0;
//
//
dz = (cc*(bbp - bb) / hp - bb*(ccp - cc) / hp) / ((bbz - bb) / hz*(ccp - cc) / hp - (bbp - bb) / hp*(ccz - cc) / hz);
dp = (cc*(bbz - bb) / hz - bb*(ccz - cc) / hz) / ((ccz - cc) / hz*(bbp - bb) / hp - (bbz - bb) / hz*(ccp - cc) / hp);
//
//
eta = eta + dz;
epsilon = epsilon + dp;
//
//
dd = fabs(bb) + fabs(cc);
//
cout << "iter = " << iter << " dd= " << dd << endl;
//
} while (dd>EPS);
//
//
//
u[0] = 0.0;
y[0] = 0.0;
z[0] = eta;
w[0] = 1.0;
p[0] = epsilon;
//
//
for (i = 0; i<n; i++)
{
a1 = f1(x[i], u[i], y[i], z[i], w[i], p[i])*h;
b1 = f2(x[i], u[i], y[i], z[i], w[i], p[i])*h;
c1 = f3(x[i], u[i], y[i], z[i], w[i], p[i])*h;
d1 = f4(x[i], u[i], y[i], z[i], w[i], p[i])*h;
g1 = f5(x[i], u[i], y[i], z[i], w[i], p[i])*h;
//
a2 = f1(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1)*h;
b2 = f2(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1)*h;
c2 = f3(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1)*h;
d2 = f4(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1)*h;
g2 = f5(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1)*h;
//
a3 = f1(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2)*h;
b3 = f2(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2)*h;
c3 = f3(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2)*h;
d3 = f4(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2)*h;
g3 = f5(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2)*h;
//
a4 = f1(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3)*h;
b4 = f2(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3)*h;
c4 = f3(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3)*h;
d4 = f4(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3)*h;
g4 = f5(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3)*h;
//
u[i + 1] = u[i] + (a1 + 2.0*a2 + 2.0*a3 + a4) / 6.0;
y[i + 1] = y[i] + (b1 + 2.0*b2 + 2.0*b3 + b4) / 6.0;
z[i + 1] = z[i] + (c1 + 2.0*c2 + 2.0*c3 + c4) / 6.0;
w[i + 1] = w[i] + (d1 + 2.0*d2 + 2.0*d3 + d4) / 6.0;
p[i + 1] = p[i] + (g1 + 2.0*g2 + 2.0*g3 + g4) / 6.0;
}
//
//
for (i = 0; i<n + 1; i++)
{
printf("x[i] = %lf then u[i]= %lf\n", x[i], u[i]);
printf("x[i] = %lf then y[i]= %lf\n", x[i], y[i]);
printf("x[i] = %lf then z[i]= %lf\n", x[i], z[i]);
printf("x[i] = %lf then w[i]= %lf\n", x[i], w[i]);
printf("x[i] = %lf then p[i]= %lf\n", x[i], p[i]);

}
delete[]x;
delete[]u;
delete[]y;
delete[]z;
delete[]w;
delete[]p;
//
system("pause");
}">

What I have tried:



I have tried to change \"eta\" and \"epsilon\" but nothing special happens.

What I have tried:

I have tried to change "eta" and "epsilon" but nothing special happens.

推荐答案

You should learn to use the debugger as soon as possible. Rather than guessing what your code is doing, It is time to see your code executing and ensuring that it does what you expect.



The debugger allow you to follow the execution line by line, inspect variables and you will see that there is a point where it stop doing what you expect.

Debugger - Wikipedia, the free encyclopedia[^]

Mastering Debugging in Visual Studio 2010 - A Beg inner’s Guide[^]

http://docs.oracle.com/javase/7/docs/technotes/tools/windows/jdb.html[^]

https://www.jetbrains.com/idea/help/debugging-your-first-java-application.html[^]



With the debugger, you will see what your code is really doing.

You should learn to use the debugger as soon as possible. Rather than guessing what your code is doing, It is time to see your code executing and ensuring that it does what you expect.

The debugger allow you to follow the execution line by line, inspect variables and you will see that there is a point where it stop doing what you expect.
Debugger - Wikipedia, the free encyclopedia[^]
Mastering Debugging in Visual Studio 2010 - A Beginner's Guide[^]
http://docs.oracle.com/javase/7/docs/technotes/tools/windows/jdb.html[^]
https://www.jetbrains.com/idea/help/debugging-your-first-java-application.html[^]

With the debugger, you will see what your code is really doing.

这篇关于用runge kutta方法求解颂歌系统的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！