求变步长四阶龙格库塔法的c程序！！
一直困扰小弟，这里就拜各位大虾了！！:handshake 请参阅《C常用算法程序集》第二版的P177；或第三版的P294.如果找不到这两本书，发短消息，或邮件konggu168@sohu.com /* vsrk4.c: Fourth-order Runge-kutta method with variable step. */#define N 100
#include "vsrk4.c"
double f( double x, double y );
extern int vsrk4( double hm, double x01,
double x, double y, double eps);
void main( )
{
int i,k;
static double y0,eps,x01,hm,x,y;
eps=1.0e-8;
x01=0.0; /* x01=a */
x01=1.0; /* x01=b */
hm=0.0005; hm=0.1;
y0=1.0; y=y0;
x=x01;
system("cls");
printf("Variable step Fourth-order Runge-Kutta method:\n");
k=vsrk4( hm, x01, x, y, eps );
printf(" xn yn\n");
for(i=0;i<=k;i++)
printf("%12.8f %12.8f\n",x,y);
getch();
}
/********* function f( x, y ) ********/
double f( double x, double y )
{
double f;
f=y-2.0*x/y;
return(f);
}
/* vsrk4.c: Fourth-order Runge-kutta method with variable step. */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
extern double f( double x, double y );
double rk4( double x, double y, double hh );
int vsrk4(double hm,double x01,double x[],double y[],double eps)
{
int i;
double xx,h,hmax,hmin,x0,x1,yn,ynh,delta;
hmin=hm; hmax=hm;
x0=x01;x1=x01;
h=hmax;
xx=x0;
i=0;
x=xx;
while( xx < x1 )
{
again:if( xx+h > x1 ) h=x1-xx;
yn=rk4(xx,y,h);
ynh=rk4(xx,y,h/2.0);
ynh=rk4(xx+h/2.0,ynh,h/2.0);
delta=fabs((yn-ynh)/15.0);
if( ( delta < eps )||( h == hmin ) )
{
xx=xx+h;
i=i+1;
x=xx;
y=ynh;
if( i >= 99 ) goto endd;
if( delta < 0.05*eps )
{
h=2.0*h;
if( h > hmax ) h=hmax;
}
}
else
{
h=h/2.0;
if( h < hmin ) h=hmin;
goto again;
}
}
endd:return(i);
}
/**************** fourth order Runge-Kutta *************************/
double rk4( double xn, double yn, double h )
{
double ynp,k1,k2,k3,k4;
k1=h*f(xn,yn);
k2=h*f(xn+0.5*h,yn+0.5*k1);
k3=h*f(xn+0.5*h,yn+0.5*k2);
k4=h*f(xn+h,yn+k3);
ynp=yn+(k1+2.0*k2+2.0*k3+k4)/6.0;
return(ynp);
} 回复 3 # 风花雪月 的帖子
大侠，有变步长龙格库塔法的matlab程序吗？急求！不胜感激 回复 4 # 煜宸0922 的帖子
Matlab还用编写RK程序吗？软件自带的ode45比起任何已经共享的C程序来说，精确多了，那可是一帮专门搞算法的大侠编写的供调用的function啊 我也来蹭个程序学学。。。 回复 3 # 风花雪月 的帖子
这个代码有误。。。
页:
[1]