/*本程序需要输入初始值r，多少都可以。但是会有精度上的差别*/
#include "stdio.h"
#include "math.h"
main()
{
float x,y,u,u1,u2,u10,u20,r,h,w,t,wf,wff,s1,s2,s3,w0,error1,error2,ht;
printf("Please input r=");
scanf("%f",&r);
h=0.00005;
ht=0.1;
for(r=r;r<=20.0;r=r+2.0)
{
for(x=0.0,u10=0.0;x<=r;x=x+4*h)
{
u1=2.0*h*(7.0*0.5*x*x*exp((-1.0)*x)+32.0*0.5*(x+h)*(x+h)*exp((-1.0)*(x+h))+12.0*0.5*(x+2.0*h)*(x+2.0*h)*exp((-1.0)*(x+2.0*h))+32.0*0.5*(x+3.0*h)*(x+3.0*h)*exp((-1.0)*(x+3.0*h))+7.0*0.5*(x+4.0*h)*(x+4.0*h)*exp((-1.0)*(x+4.0*h)))/45.0;
u10=u10+u1;
}
for(y=0.000000000000000001,u20=0.0;y<=1/r+0.000000000000000001;y=y+4*h)
{
u2=2.0*h*(7.0*0.5*(1.0/y)*(1.0/y)*(1.0/y)*exp((-1.0)/y)+32.0*0.5*(1.0/(y+h))*(1.0/(y+h))*(1.0/(y+h))*exp((-1.0)/(y+h))+12.0*0.5*(1.0/(y+2.0*h))*(1.0/(y+2.0*h))*(1.0/(y+2.0*h))*exp((-1.0)/(y+2.0*h))+32.0*0.5*(1.0/(y+3.0*h))*(1.0/(y+3.0*h))*(1.0/(y+3.0*h))*exp((-1.0)/(y+3.0*h))+7.0*0.5*(1.0/(y+4.0*h))*(1.0/(y+4.0*h))*(1.0/(y+4.0*h))*exp((-1.0)/(y+4.0*h)))/45.0;
u20=u20+u2;
}
u=u10+r*u20;
for(t=0.0,wff=0,wf=1-0.5*(ht+2.0)*exp((-1.0)*ht);t<=r-2.0*ht;t=t+ht)
{
s1=(-0.5)*t*exp((-1.0)*t);
s2=(-0.5)*(t+ht)*exp((-1.0)*(t+ht));
s3=(-0.5)*(t+2.0*ht)*exp((-1.0)*(t+2.0*ht));
w=2.0*wf-wff+ht*ht*(s3+10.0*s2+s1)/12.0;
wff=wf;
wf=w;
}
w0=1.0-0.5*(r+2.0)*exp((-1.0)*r);
error1=w0-u;
error2=w0-w;
printf("exact W=%f ",w0);
printf("green function error=%f ",error1);
printf("directly error=%f\n",error2);
}
}
