自己编写4阶Runge-Kutta法的时候，调试有问题，麻烦大家看看！

jungwoo

里面我用了一个非常简单的函数f(t,u)=2*t+u,初值为u（0）=1，范围为（0，0.5）
代码：
// RK4.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include "stdio.h"
#include "stdlib.h"
#include "math.h"


class RK
{
        private:
                double k1,k2,k3,k4;
                double h,b,a,u;
                double *result_u;
                FILE * fp1;
        public:
                void seth(double l=0)
                {
                        h=l;
                } //设步长
                void setf(double xa=0,double xb=0,double y=0) //设初值和范围（xa,xb）
                {
                        b=xb;
                        a=xa;
                        u=y;
                }
                double f(double t,double u) //函数值，修改它以适应各自需要
                {
                        //函数设定
                        double f=2*t+u;
                        return f;
                }
                void dork() //R-K 主函数
                {
                        result_u=new double[(b-a)/h+1];
                        result_u[0]=u;
                        int count;
                        for(count=0;count<(b-a)/h;count++)
                        {
                                k1=h*f(a+count*h,u);
                                k2=h*f(a+count*h+h/2,u+k1/2);
                                k3=h*f(a+count*h+h/2,u+h*k2/2);
                                k4=h*f(a+count*h+h,u+k3);
                                u=u+(k1+2*k2+2*k3+k4)/6;
                                result_u[count+1]=u;
                                printf("%d\t",u);
                        }
                        if((fp1=fopen("result_u.txt","w"))==NULL)                //如果未打开则显示Cannot open file。
                        {
                                printf("Cannot open file. \n");
                                exit(1);
                        }
                        printf("result_u:\n");
                        for(count=0;count<(b-a)/h+1;count++)
                        {
                                printf("%d\t",result_u[count]);
                                fprintf(fp1,"%d\t",result_u[count]);
                        }
                        fclose (fp1);
                }
};

void main()
{
        RK my;
        my.seth(0.1);
        my.setf(0,0.5,1.0);
        my.dork();
}

[ 本帖最后由 jungwoo 于 2008-4-14 21:25 编辑 ]

jungwoo

有人看看吗？错误找不出来，调试的时候出现的结果很奇怪