定积分

辛普生求积分：

#include <stdio.h>

#include <math.h>

 

doubleFunction(double);

doubleSIMP1(double,double,int);

doubleSIMP2(double,double,double);

 

voidmain()

{

         doublea1,b1,eps;

         intn1;

         doubleResult1;

         doubleResult2;

         a1 = 0.0;

         b1 = 0.8;

         n1 = 4;

         eps = 5e-7;

         Result1 = SIMP1(a1,b1,n1);

         Result2 = SIMP2(a1,b1,eps);

         printf("利用定步长辛普生积分结果为：/n");

         printf("I1 = %.10f/n",Result1);

         printf("利用变步长辛普生积分结果为：/n");

         printf("I2 = %e/n",Result2);

}

 

doubleSIMP1(a,b,n)

doublea;

doubleb;

intn;

{

         inti;

         doubleh,s;

         h=(a-b)/(2*n);

         s=0.5*(Function(a)-Function(b));

         for(i=1;i<=n;i++)

                  s+=2*Function(a+(2*i-1)*h)+Function(a+2*i*h);

         return((b-a)*s/(3*n));

}

 

doubleSIMP2(a,b,eps)

doublea;

doubleb;

doubleeps;

{

         intk,n;

         doubleh,t1,t2,s1,s2,p,x;

         n=1;

         h=b-a;

         t1=h*(Function(a)+Function(b))/2;

         s1 = t1;

         while(1)

         {

                   p = 0;

                  for(k=0;k<=n;k++)

                   {

                            x = a+(k+0.5)*h;

                            p+=Function(x);

                   }

                  t2=(t1+h*p)/2;

                  s2=(4*t2-t1)/3;

                  if(fabs(s2-s1)>=eps)

                   {

                            t1=t2;

                            n=n+n;

                            h=h/2;

                            s1=s2;

                            continue;

                   }

                   break;

         }

         return(s2);

}

 

doubleFunction(doublex)

{

         return(cos(x));

}

欧拉方程法：

 

#include "stdio.h"

#include "stdlib.h"

#include <math.h>

 

intFunc(y,d)

doubley[],d[];

{

         d[0]=y[1];          /*y0'=y1*/

         d[1]=-y[0];                  /*y1'=y0*/

         d[2]=-y[2];                  /*y2'=y2*/

         return(1);

}

 

voidEuler1(t,y,n,h,k,z)

intn;                   /*整型变量，微分方程组中方程的个数，也是未知函数的个数*/

intk;                   /*整型变量。积分步数(包括起始点这一步)*/

doublet;            /*双精度实型变量,对微分方程进行积分的起始点t0*/

doubleh;           /*双精度实型变量。积分步长*/

doubley[];        /*双精度实型一维数组，长度为n。存放n个未知函数yi在起始点t0处的函数值*/

doublez[];         /*双精度实型二维数组，体积为nxk。返回k个积分点(包括起始点)上的未知函数值*/

{

         externintFunc();

         inti,j;

         double *d;

         d=malloc(n*sizeof(double));

         if(d == NULL)

         {

                   printf("内存分配失败/n");

                   exit(1);

         }

         /*将方程组的初值赋给数组z[i*k]*/

         for (i=0; i<=n-1; i++)

                  z[i*k]=y[i];

         for (j=1; j<=k-1; j++)

         {

                  Func(y,d);                            /*求出f(x)*/

                   for(i=0; i<=n-1; i++)

                            y[i]=z[i*k+j-1]+h*d[i];

                  

                  Func(y,d);

                 for (i=0; i<=n-1; i++)

                            d[i]=z[i*k+j-1]+h*d[i];

                 for (i=0; i<=n-1; i++)

                   {

                            y[i]=(y[i]+d[i])/2.0;

                            z[i*k+j]=y[i];

                   }

         }

    free(d);

         return;

}

voidEuler2(t,h,y,n,eps)

intn;

doublet,h,eps,y[];

{

         inti,j,m;

         doublehh,p,q,*a,*b,*c,*d;

    a=malloc(n*sizeof(double));

    b=malloc(n*sizeof(double));

    c=malloc(n*sizeof(double));

    d=malloc(n*sizeof(double));

    hh=h;

         m=1;

         p=1.0+eps;

         for (i=0; i<=n-1; i++) a[i]=y[i];

    while (p>=eps)

    {

                 for (i=0; i<=n-1; i++)

        {

                            b[i]=y[i];

                            y[i]=a[i];

                   }

        for (j=0; j<=m-1; j++)

        {

                            for (i=0; i<=n-1; i++)

                                     c[i]=y[i];

            Func(y,d);

            for (i=0; i<=n-1; i++)

                y[i]=c[i]+hh*d[i];

            Func(y,d);

            for (i=0; i<=n-1; i++)

                d[i]=c[i]+hh*d[i];

            for (i=0; i<=n-1; i++)

                y[i]=(y[i]+d[i])/2.0;

        }

        p=0.0;

        for (i=0; i<=n-1; i++)

        {

                            q=fabs(y[i]-b[i]);

            if (q>p) p=q;

        }

        hh=hh/2.0; m=m+m;

         }

    free(a);

         free(b);

         free(c);

         free(d);

         return;

}

main()

{

         inti,j;

         doubley[3],z[3][11],t,h,x,eps;

         y[0]=-1.0;                    /*初值y0(0)=-1.0*/

         y[1]=0.0;                      /*初值y1(0)=-1.0*/

         y[2]=1.0;                      /*初值y2(0)=-1.0*/

         t=0.0;                                     /*起始点t=0*/

         h=0.01;                                  /*步长为0.01*/

         eps = 0.00001;

         Euler1(t,y,3,h,11,z);

         printf("定步长欧拉法结果：/n");

         for (i=0; i<=10; i++)

         {

                   x=i*h;

                  printf("t=%5.2f/t   ",x);

                   for(j=0; j<=2; j++)

                            printf("y(%d)=%e ",j,z[j][i]);

                  printf("/n");

         }

         y[0]=-1.0;                    /*重新赋初值*/

         y[1]=0.0;                     

         y[2]=1.0;                     

         printf("变步长欧拉法结果：/n");

         printf("t=%5.2f/t   ",t);

        for (i=0; i<=2; i++)

                  printf("y(%d)=%e ",i,y[i]);

         printf("/n");

         for (j=1; j<=10; j++)

         {

                  Euler2(t,h,y,3,eps);

                  t=t+h;

                  printf("t=%5.2f/t   ",t);

                 for (i=0; i<=2; i++)

                            printf("y(%d)=%e ",i,y[i]);

                  printf("/n");

         }

}