//IEulerRK.cpp--Improved Euler and Runge-Kutta(4)
//qiu changweifen fangcheng shuzhijie
#include < stdio.h >
#include < stdlib.h >
#include < math.h >
#define FMT "%-9.5g"
typedef double dbl;
//prototypesdbl
fxy(dbl x, dbl y);
dbl f(dbl x);
void IEuler(dbl x0, dbl y0, dbl b, dbl h);
void RK(dbl x0, dbl y0, dbl b, dbl h);
int main() {
printf("\nImproved Euler:\n");
IEuler(0.0, 1.0, 1.0, 0.1);
printf("\nRunge-Kutta:\n");
RK(0.0, 1.0, 1.0, 0.2);
return 0;
}
dbl fxy(dbl x, dbl y) {
if (fabs(y) > 1e-7) return (y - 2 * x / y);
else exit( - 2);
return 0.0;
}
dbl f(dbl x) {
if (1 + 2 * x >= 0.0) return sqrt(1 + 2 * x);
else exit( - 2);
return 0.0;
}
void IEuler(dbl x0, dbl y0, dbl b, dbl h) {
dbl yp,
yc,
xn,
yn;
printf("xn yn y(xn)\n");
xn = x0;
yn = y0;
while (xn <= b - 1e-7) {
yp = yn + h * fxy(xn, yn);
yc = yn + h * fxy(xn + h, yp);
yn = (yp + yc) / 2.0;
xn += h;
printf(FMT, xn);
printf(FMT, yn);
printf(FMT, f(xn));
printf("\n");
}
}
void RK(dbl x0, dbl y0, dbl b, dbl h) {
dbl k1,
k2,
k3,
k4;
dbl xn,
yn;
printf("xn yn y(xn)\n");
xn = x0;
yn = y0;
while (xn <= b - 1e-7) {
k1 = fxy(xn, yn);
k2 = fxy(xn + h / 2.0, yn + k1 * h / 2.0);
k3 = fxy(xn + h / 2.0, yn + k2 * h / 2.0);
k4 = fxy(xn + h, yn + k3 * h);
yn = yn + (k1 + 2 * k2 + 2 * k3 + k4) * h / 6;
xn += h;
printf(FMT, xn);
printf(FMT, yn);
printf(FMT, f(xn));
printf("\n");
}
}
IEulerRK
原创
©著作权归作者所有:来自51CTO博客作者mb643d15e043b20的原创作品,请联系作者获取转载授权,否则将追究法律责任
提问和评论都可以,用心的回复会被更多人看到
评论
发布评论
相关文章