in this code have a few of function .Function of lagrange, euler, interpolate delay, and the function to get the exact solution..
#include <stdio.h>#include <math.h>#define PI 4.0E0*atan(1.0E0)#define IP 1000#define TRUE 1#define FALSE 0double XOUT[IP],YOUT[IP],ZOUT[IP];double f[IP],X[IP],Yx[IP],Yalpha[IP];double FY(),FALPHA(),FF(),FSOLN(),FYALPHA(),FINTERPOLATE1(),FINTERPOLATE2();void FEULER();main() { double H,ENDPOINT,EC,alpha,Ypredict,Ycorrect,soln,error; int I,NEXT; H = 0.01; X[0] = 0.0;//initial value ENDPOINT=1.00; EC = 0.000001; NEXT=TRUE; //compute initial values Yx[0]=FY(X[0]); alpha=FALPHA(X[0],Yx[0]); Yalpha[0]=FY(alpha); f[0]=FF(X[0],Yx[0],Yalpha[0]); //store initial values XOUT[0]=X[0]; YOUT[0]=Yx[0]; ZOUT[0]=f[0]; FEULER(H); I=3; X[I+1]=X[I]+H; printf("I=%d\n",I); while(NEXT){ Ypredict=Yx[I]+H/24.0*(55.0*f[I]-59.0*f[I-1]+37.0*f[I-2]-9*f[I-3]); Yalpha[I+1]=FYALPHA(I,H); p1:f[I+1]=FF(X[I+1],Ypredict,Yalpha[I+1]); Ycorrect=Yx[I]+H/24.0*(9.0*f[I+1]+19.0*f[I]-5.0*f[I-1]+f[I-2]); if (fabs(Ypredict - Ycorrect) > EC){ Ypredict=Ycorrect; goto p1; } f[I+1]=FF(X[I+1],Ycorrect,Yalpha[I+1]); Yx[I+1]=Ycorrect; X[I+1]=X[I]+H; //store values XOUT[I+1]=X[I+1]; YOUT[I+1]=Yx[I+1]; ZOUT[I+1]=f[I+1]; //error calculation soln=FSOLN(X[0]); error=fabs(Yx[I+1]-soln); printf("X[%d]=%9.6f Y=%9.6e Ytrue=%9.6e error=%9.6e\n", I+1,X[I+1],Yx[I+1],soln,error); //check if endpoint is reached if(X[I+1]>ENDPOINT) NEXT=FALSE; else I=I+1; }//end while}//end main//initial functiondouble FY(T)double T;{ double Y; Y=sin(T); return Y;}//function alphadouble FALPHA(T,YT)double T,YT;{ double ALPHA; ALPHA=T-PI; return ALPHA;}//function Fdouble FF(T,YT,YA)double T,YT,YA;{ double F; F=-YA; return F;}//function exact solutiondouble FSOLN(T)double T;{ double SOLN; SOLN=sin(T); return SOLN;}//function yalphadouble FYALPHA(I,H)int I;double H;{ int J,N; double ALP,YALPHA,S; ALP=FALPHA(X[I+1],Yx[I+1]); if(ALP<=X[0]) YALPHA=FY(ALP); else if((ALP>X[3])&&(ALP<=X[4])){ YALPHA=FINTERPOLATE2(I,ALP); } else if((ALP>X[0])&&(ALP<=X[3])){ J=0; while ((ALP<=X[J])||(ALP>X[J+1])) J=J+1; S=(ALP-X[J])/H; YALPHA=YOUT[J]+S*H*ZOUT[J]; } else{ J=3; while((ALP<=X[J])||(ALP>X[J+1])) J=J+1; N=5; YALPHA=FINTERPOLATE1(J,N,ALP); } return YALPHA;}//function interpolate delay1double FINTERPOLATE1(J,N,ALP)int J,N;double ALP;{ int d,i,j; double x[15],y[15],D[15],Q[15][15]; double DELAY; d=N-2; for(i=0;i<=N;i++){ x[i]=XOUT[J-d+i]; y[i]=YOUT[J-d+i]; Q[i][0]=y[i]; } D[0]=ALP-x[0]; for(i=1;i<=N;i++){ D[i]=ALP-x[i]; for(j=1;j<=i;j++) Q[i][j]=(D[i]*Q[i-1][j-1]-D[i-j]*Q[i][j-1])/(D[i]-D[i-j]); } DELAY=Q[N][N]; return DELAY;}//function interpolate delay1double FINTERPOLATE2(N,ALP)int N;double ALP;{ int i,j,NP; double x[15],y[15],D[15],Q[15][15]; double DELAY; if(N<4) NP=N; else NP=4; for(i=0;i<=NP;i++){ x[i]=XOUT[N-NP+i]; y[i]=YOUT[N-NP+i]; Q[i][0]=y[i]; } D[0]=ALP-x[0]; for(i=1;i<=NP;i++){ D[i]=ALP-x[i]; for(j=1;j<=i;j++) Q[i][j]=(D[i]*Q[i-1][j-1]-D[i-j]*Q[i][j-1])/(D[i]-D[i-j]); } DELAY=Q[NP][NP]; return DELAY;}//function EULERvoid FEULER(H)double H;{ int p,j; double alpha,r; for(p=1;p<=3;p++){ X[p]=X[0]+p*H; //store X XOUT[p]=X[p]; } for(p=1;p<=3;p++){ Yx[p]=Yx[p-1]+H*f[p-1]; alpha=FALPHA(X[p],Yx[p]); if(alpha<=X[0]) Yalpha[p]=FY(alpha); else if((alpha>X[p-1])&&(alpha<=X[p])){ r=(alpha-X[p-1])/H; Yalpha[p]=Yx[p-1]+r*H*f[p-1]; } else{ j=0; while((alpha<=X[j])||(alpha>X[j+1])) j=j+1; r=(alpha-X[j])/H; Yalpha[p]=Yx[j-1]+r*H*f[j-1]; } f[p]=FF(X[p],Yx[p],Yalpha[p]); //store computed values YOUT[p]=Yx[p]; ZOUT[p]=f[p]; }}
Best Answer