标题:学过计算方法的帮忙看看吧,发出来也没人看一下
只看楼主
梁朝斌
Rank: 4
等 级:业余侠客
帖 子:192
专家分:288
注 册:2012-10-21
结帖率:100%
 问题点数:0 回复次数:0 
学过计算方法的帮忙看看吧,发出来也没人看一下

把不同的方法合起来写,核心算法完全正确,请各位帮忙,我那里不对,还要怎么改
 下面这几个方法都是单独验证的时候正确的,也能得到完全正确的解,步长step=0.001,区间a=0,b=0.01,核心算法就是类跟:dy/dx=x+y;y=1+exp(x)-x,通过微分求得y最终的解,子程序中f1是微分方程,f2是准确的解就是这么个意思,但是一天了也没人能给说说
#include<stdio.h>
#include<math.h>
#define maxsize 100
double f1(double x,double y);
double f2(double x);
/*
** 欧拉折线法求解微分方程
**/
void eulerian_method(double x[maxsize],double y[maxsize])
{
    double a,b,step;
    double f1(double x,double y);
    double f2(double x);
    int i,n;
    printf("please input step=\n");/*输入步长*/
    scanf("%lf",&step);

    printf("please input a= and b=\n");/*微分区间*/
    scanf("%lf %lf",&a,&b);
    n=(int)((b-a)/step);

    for(i=0;i<=n;i++)/*分区间*/
        x[i]=a+i*step;

    y[0]=0;
    for(i=0;i<=n-1;i++)
    {
     y[i+1]=y[i]+step*f1(x[i+1],y[i]);/*算法核心,求解公式*/
    }
    for(i=0;i<=n;i++)
    {
        printf("x[%d]=%lf,y[%d]=%lf,f2=%lf",i,x[i],i,y[i],f2(x[i]));/*输出x[i],y[i],准确值*/
        printf("\n");
    }
}
/*
** 欧拉折线法改进求解微分方程
**/
void modified_eulerian_method(double x[maxsize],double y[maxsize])
{
    double a,b,step,k1,k2;
    double f1(double x,double y);
    double f2(double x);
    int i,n;

    printf("please input step=\n");/*输入步长*/
    scanf("%lf",&step);

    printf("please input a= and b=\n");/*微分区间*/
    scanf("%lf %lf",&a,&b);

    n=(int)((b-a)/step);/*分区间*/
    for(i=0;i<=n;i++)
        x[i]=a+i*step;
    y[0]=0;
    for(i=0;i<=n-1;i++)
    {
        k1=step*f1(x[i],y[i]);
        k2=step*f1(x[i]+step,y[i]+k1);/*算法核心,求解公式*/
        y[i+1]=y[i]+(k1+k2)/2;
    }
    for(i=0;i<=n;i++)
    {
        printf("x[%d]=%lf,y[%d]=%lf,f2=%lf",i,x[i],i,y[i],f2(x[i]));/*输出x[i],y[i],准确值*/
        printf("\n");
    }
}
/*
** 龙格-库塔公式求解微分方程
**/
void runge_kutta(double x[maxsize],double y[maxsize])
{
    double a,b,step,k1,k2,k3,k4;
    double f1(double x,double y);
    double f2(double x);
    int i,n;
    printf("please input step=\n");/*输入步长*/
    scanf("%lf",&step);

    printf("please input a= and b=\n");/*微分区间*/
    scanf("%lf %lf",&a,&b);

    n=(int)((b-a)/step);
    for(i=0;i<=n;i++)
        x[i]=a+i*step;
    y[0]=0;
    for(i=0;i<=n-1;i++)
    {
           k1=step*f1(x[i],y[i]);
        k2=step*f1(x[i]+step/2,y[i]+k1/2);
        k3=step*f1(x[i]+step/2,y[i]+k2/2);/*算法核心,求解公式*/
        k4=step*f1(x[i]+step,y[i]+k3);
        y[i+1]=y[i]+(k1+2*k2+2*k3+k4)/6;
    }
    for(i=0;i<=n;i++)
    {
        printf("x[%d]=%lf,y[%d]=%lf,f2=%lf",i,x[i],i,y[i],f2(x[i]));
        printf("\n");
    }
}
/*
** 基尔公式-求解微分方程
**/
void gill(double x[maxsize],double y[maxsize])
{
    double a,b,step,k1,k2,k3,k4;
    double f1(double x,double y);
    double f2(double x);
    int i,n;
    printf("please input step=\n");/*输入步长*/
    scanf("%lf",&step);

    printf("please input a= and b=\n");/*微分区间*/
    scanf("%lf %lf",&a,&b);

    n=(int)((b-a)/step);
    for(i=0;i<=n;i++)
        x[i]=a+i*step;
    y[0]=0;
    for(i=0;i<=n-1;i++)
    {
           k1=f1(x[i],y[i]);
        k2=f1(x[i]+step/2,y[i]+k1*step/2);
        k3=f1(x[i]+step/2,y[i]+(sqrt(2)-1)*step*k1/2+(1-sqrt(2)/2)*step*k2);/*算法核心,求解公式*/
        k4=f1(x[i]+step,y[i]-sqrt(2)*step*k2/2+(1+sqrt(2)/2)*step*k3);
        y[i+1]=y[i]+(k1+(2-sqrt(2))*k2+(2+sqrt(2))*k3+k4)*step/6;
    }
    for(i=0;i<=n;i++)
    {
        printf("x[%d]=%lf,y[%d]=%lf,f2=%lf",i,x[i],i,y[i],f2(x[i]));/*输出x[i],y[i],准确值*/
        printf("\n");
    }
}
/*下面两个是子程序*/
double f1(double x,double y)
{
    double z;
    z=622*sin(314*x)-20*y;
    return z;
}
double f2(double x)
{
    double z;
    z=48827*(exp(-20*x)+10*sin(314*x)/157-cos(314*x))/24749;
    return z;
}

void menu()
{
    printf("********************************\n");
    printf("请选择命令号!                  *\n");
    printf("1欧拉折线法求解微分!           *\n");
    printf("2欧拉改进法求解微分!           *\n");
    printf("3龙格-库塔公式求解微分!        *\n");
    printf("4-基尔公式求解微分!            *\n");
    printf("0退出!                         *\n");
    printf("********************************\n");
}
int main(void)
{
    int sel=0;
    double x[maxsize],y[maxsize];
    while(sel!=0)
    {
        menu();
        printf("请输入命令号:\n");
        scanf("%d",&sel);
        switch(sel)
        {
            case 1:eulerian_method(x,y);
                break;

            case 2:modified_eulerian_method(x,y);
                break;

            case 3: runge_kutta(x,y);
                break;

            case 4:gill(x,y);
                break;

            case 0:exit(1);
                break;

            default:
                 printf("输入的命令号错误!请重新输入:\n");
                 break;
        }
    }
    return 0;
}

搜索更多相关主题的帖子: double include 计算方法 子程序 
2012-12-16 17:03



参与讨论请移步原网站贴子:https://bbs.bccn.net/thread-391503-1-1.html




关于我们 | 广告合作 | 编程中国 | 清除Cookies | TOP | 手机版

编程中国 版权所有,并保留所有权利。
Powered by Discuz, Processed in 0.073115 second(s), 8 queries.
Copyright©2004-2024, BCCN.NET, All Rights Reserved