积分和微分实现方法
#include#include
#include
#define eps 0.000001
#define maxsize 100
#pragma warning(disable:4996)
/*
**f1是微分公式
*/
double f1(double x,double y)
{
double z;
z=622*sin(314*x)-20*y;
return z;
}
/*
**f2是准确值的公式
*/
double f2(double x)
{
double z;
z=48827*(exp(-20*x)+10*sin(314*x)/157-cos(314*x))/24749;
return z;
}
/*
**f3是要求积分的公式
*/
double f3(double x)
{
double z;
z=sqrt(1+pow(x,2));
return z;
}
/*
**欧拉折线法求y的值
*/
void eulerian_method(void)
{
double a,b,step;
int i,n;
double x[maxsize], y[maxsize];
printf("请输入步长step=");
scanf("%lf",&step);
printf("请输入积分下限a=");
scanf("%lf",&a);
printf("请输入积分上限b=");
scanf("%lf",&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]));
printf("\n");
}
}
/*
**改进欧拉法求y的值
*/
void modified_eulerian_method(void)
{
double a,b,step,k1,k2;
int i,n;
double x[maxsize], y[maxsize];
printf("请输入步长step=");
scanf("%lf",&step);
printf("请输入积分下限a=");
scanf("%lf",&a);
printf("请输入积分上限b=");
scanf("%lf",&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]));
printf("\n");
}
}
/*
**龙格-库塔公式求y的值
*/
void runge_kutta(void)
{
double a,b,step,k1,k2,k3,k4;
int i,n;
double x[maxsize],y[maxsize];
printf("请输入步长step=");
scanf("%lf",&step);
printf("请输入积分下限a=");
scanf("%lf",&a);
printf("请输入积分上限b=");
scanf("%lf",&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");
}
}
/*
**基尔公式求y的值
*/
void gill(void)
{
double a,b,step,k1,k2,k3,k4;
int i,n;
double x[maxsize],y[maxsize];
printf("请输入步长step=");
scanf("%lf",&step);
printf("请输入积分下限a=");
scanf("%lf",&a);
printf("请输入积分上限b=");
scanf("%lf",&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.0)-1)*step*k1/2+(1-sqrt(2.0)/2)*step*k2);
k4=f1(x[i]+step,y[i]-sqrt(2.0)*step*k2/2+(1+sqrt(2.0)/2)*step*k3);
y[i+1]=y[i]+(k1+(2-sqrt(2.0))*k2+(2+sqrt(2.0))*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]));
printf("\n");
}
}
/*
**梯形公式求y的积分
*/
void trapezoid_formula(void)
{
double a,b,t,x[maxsize],step,s0,s1;
double f3(double x);
int i,n;
printf("请输入积分下限a=");
scanf("%lf",&a);
printf("请输入积分上限b=");
scanf("%lf",&b);
n=1;
s0=(f3(a)+f3(b))*(b-a)/2;
do
{
n=n+1;
step=(b-a)/n;
for(i=0;i<=n;i++)
x[i]=a+i*step;
s1=0;
for(i=0;i<=n-1;i++)
s1+=(f3(x[i])+f3(x[i+1]))*step/2;
t=s0;
s0=s1;
}
while (fabs(s0-t)>eps);
{
printf("求得积分s0=%lf\n",s0);
printf("此时步长为:%d",n);
printf("\n");
}
}
/*
**辛普生公式求y的积分
*/
void simpson(void)
{
double a,b,t,x[maxsize],step,s0,s1;
double f3(double x);
int i,n;
printf("请输入积分下限a=");
scanf("%lf",&a);
printf("请输入积分上限b=");
scanf("%lf",&b);
n=0;
s0=(f3(a)+f3(b)+f3((a+b)/2))*(b-a)/6;
do
{
n=n+1;
step=(b-a)/n;
for(i=0;i<=n;i++)
x[i]=a+i*step;
s1=0;
for(i=0;i<=n-1;i++)
s1+=(f3(x[i])+f3(x[i+1])+4*f3((x[i]+x[i+1])/2))*step/6;
t=s0;
s0=s1;
}
while (fabs(s0-t)>eps&&n<100);
{
printf("求得积分s0=%lf\n",s0);
printf("此时步长为:%d",n);
printf("\n");
}
}
/*
**科特斯公式求y的积分
*/
void Cortez_formula(void)
{
double a,b,t,x[maxsize],step,s0,s1;
double m1,m2,m3,sum1,sum2,sum3,sum4;
double f3(double x);
int i,n;
printf("请输入积分下限a=");
scanf("%lf",&a);
printf("请输入积分上限b=");
scanf("%lf",&b);
n=0;
s0=0;
do
{
n=n+1;
step=(b-a)/n;
for(i=0;i<=n;i++)
x[i]=a+i*step;
sum1=sum2=sum3=sum4=0;
for(i=0;i<=n-1;i++)
{
m1=x[i]+(x[i+1]-x[i])/4;
sum1+=f3(m1);
}
for(i=0;i<=n-1;i++)
{
m2=x[i]+(x[i+1]-x[i])/2;
sum2+=f3(m2);
}
for(i=0;i<=n-1;i++)
{
m3=x[i]+(x[i+1]-x[i])*3/4;
sum3+=f3(m3);
}
for(i=0;i<=n-1;i++)
{
sum4+=f3(x[i]);
}
s1=0;
s1+=(7*(f3(a)+f3(b))+32*sum1+12*sum2+32*sum3+14*sum4)*step/90;
t=s0;
s0=s1;
}
while (fabs(s0-t)>eps&&n<100);
{
printf("求得积分s0=%lf\n",s0);
printf("\n");
}
}
/*
**龙贝格公式求y的积分
*/
void Dragon_berger(void)
{
double a,b,t[maxsize],s[maxsize],c[maxsize],r[maxsize];
double f3(double x);
double sum;
int i,n,k;
printf("请输入积分下限a=");
scanf("%lf",&a);
printf("请输入积分上限b=");
scanf("%lf",&b);
printf("please input n=\n");
scanf("%d",&n);
t[0]=(b-a)*(f3(a)+f3(b))/2;
printf("t[0]=%lf\n",t[0]);
for(k=1;k<=n-1;k++)
{
sum=0;
for(i=1;i<=pow(2,k-1);i++)
sum+=f3(a+(b-a)*(2*i-1)/pow(2,k));
t[k]=t[k-1]/2+(b-a)*sum/pow(2,k);
printf("t[%d]=%lf",k,t[k]);
printf("\n");
}
printf("\n");
for(k=0;k<=n-2;k++)
{
s[k]=4*t[k+1]/3-t[k]/3;
printf("s[%d]=%lf",k,s[k]);
printf("\n");
}
printf("\n");
for(k=0;k<=n-3;k++)
{
c[k]=16*s[k+1]/15-s[k]/15;
printf("c[%d]=%lf",k,c[k]);
printf("\n");
}
printf("\n");
for(k=0;k<=n-4;k++)
{
r[k]=64*c[k+1]/63-c[k]/63;
printf("r[%d]=%lf",k,r[k]);
printf("\n");
}
}
void menu()
{
printf("*********************************\n");
printf("*请选择命令号! *\n");
printf("*1 欧拉折线法求解微分! *\n");
printf("*2 欧拉改进法求解微分! *\n");
printf("*3 龙格-库塔公式求解微分! *\n");
printf("*4 基尔公式求解微分! *\n");
printf("*********************************\n");
printf("*5 梯形公式求积分! *\n");
printf("*6 辛普生公式求积分! *\n");
printf("*7 科特斯公式求积分! *\n");
printf("*8 龙贝格公式求积分! *\n");
printf("*0 退出! *\n");
printf("*********************************\n");
}
int main(void)
{
int sel = 0;
while(1)
{
menu();
printf("请输入命令号:\n");
scanf("%d",&sel);
if ((sel<0)||(sel>99))
{
printf("输入的命令号错误,请重新输入! n");
return 0;
}
switch(sel)
{
case 1:eulerian_method();
break;
case 2:modified_eulerian_method();
break;
case 3:runge_kutta();
break;
case 4:gill();
break;
case 5:trapezoid_formula();
break;
case 6:simpson();
break;
case 7:Cortez_formula();
break;
case 8:Dragon_berger();
break;
case 0:exit(1);
break;
default:printf("输入的命令号错误,请重新输入!");break;
}
}
return 0;
}
