Hi..I would like to convert from c to matlab..anybody can help me?? below, is the C coding.

Question

0 votes

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 0
double 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 function
double FY(T)
double T;
{
  double Y;
  Y=sin(T);
  return Y;
}
//function alpha
double FALPHA(T,YT)
double T,YT;
{
  double ALPHA;
  ALPHA=T-PI;
  return ALPHA;
}
//function F
double FF(T,YT,YA)
double T,YT,YA;
{
  double F;
  F=-YA;
  return F;
}
//function exact solution
double FSOLN(T)
double T;
{ 
  double SOLN;
   
  SOLN=sin(T);
  return SOLN;
}
//function yalpha
double 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 delay1
double 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 delay1
double 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 EULER
void 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];
  }
}