🚀 go-pugleaf

Path: inetnews.pcvan.or.jp
From: DWF78073@pcvan.or.jp (Kazutake IMANI)
Newsgroups: fj.sci.math
Subject: Milne Simutaneous Diff.Eq. Solver 
Date: 4 November, 1995 
Distribution: sci.math

$B$O$8$a$^$7$F!"%$%^%K!!:#?N!!OBIp!!(B@ $B<+M36H$G$9!#(B

$B!!(BMilne $BK!$GO"N)(B(9$B85(B)$B>oHyJ,J}Dx<0$r!"#C8@8l$G2r$/%W%m%0%i%`$r(B
$B65$($F$/$@$5$$!#(B  $B!!0l85$N>l9g$N%j%9%H$O!"2<5-$N$h$&$G$9!#(B
/* Milne method for Differential Eq.   milne.c  */
#include <stdio.h>
#include <math.h>
#include <graph.h>
#include <stdlib.h>
double rkm(double(*)(double,double),double,double,double);
double aml(double(*)(double,double),double*,double,double*,double,double);
double func(double,double);
void main ()
{static double eps,h,x,x0,xn,y[5],df[5]; static int i,j,n;
 printf("\nx0="); scanf("%le",&x0); printf("\ny="); scanf("%le", y);
 printf("\nxn="); scanf("%le",&xn); printf("\nn= "); scanf("%d",&n);
 eps=0.000001; h=(xn-x0)/n;
 _setvideomode(_98RESSCOLOR);  /* Runge Kutta  */
 for(i=0;i<=3;i++) { x=x0+i*h; df[i]=func(x,y[i]);
  printf(" %25.8e%22.8e \n",x,y[i]);  _setpixel(i,100-y[i]*50);
  if (i<3)  y[i+1]=rkm(func,x,y[i],h); }
 /*  Milne */
 for (i=4;i<=n;i++) { x=x0+i*h;
     y[4]= aml(func,df,x,y,h,eps);
     df[4]=func(x,y[4]);
     printf(" %25.8e,%22.8e,%22.8e \n",x,y[4],sin(x));
     _setpixel(i,100-50*y[4]);
     for (j=0;j<=3;j++)
        { y[j]= y[j+1]; df[j]=df[j+1]; }
 }
 exit(0);
}
/*  Runge Kutta */
double rkm(f,x,y,h)
  double (*f)(double,double);
  double x, y,h ;
{double ck0,ck1,ck2,ck3;
  ck0=(*f)(x,y)*h; ck1=(*f)(x+h*0.5,y+ck0*0.5)*h;
  ck2=(*f)(x+h*0.5,y+ck1*0.5)*h; ck3=(*f)(x+h,y+ck2)*h;
  return(y+(ck0+(ck1+ck2)*2.0+ck3)/6.0);
}
/* Milne method */
double aml(f,df,x,y,h,eps)
      double (*f)(double,double);
      double df[ ],x,y[ ],h,eps ;
{double amln,ff,py;
 py=y[0]+((df[3]+df[1])*2.0 -df[2])*h*4.0/3.0;
 while (1) {
    ff = (*f) (x,py);
    amln=y[2]+(ff+df[3]*4.0+df[2])*h/3.0;
    if (fabs(py-amln) <=eps)
      return(amln);
    py = amln;
 }
}
/*  Function */
double func (x,y)
  double x,y;
{ return( cos(x) );}

$B3'MM$+$i$N%3%a%s%H$r!"EEOC$d(B E-Mail $B$G$*BT$A$7$F$$$^$9!#(B 
$B!!(B312 $B0q>k8)$R$?$A$J$+;TElBgEg#4!]#1#6!]#4(B       TEL 029 273 3953
      PCVAN DWF78073                      $B:#?N(B $BOBIp(B $B%$%^%K%+%:(B
(B