#include "stdio.h"
#include "stdlib.h"
#include <math.h>
int Func(y,d)
double y[],d[];
{
d[0]=y[1]; /*y0'=y1*/
d[1]=-y[0]; /*y1'=y0*/
d[2]=-y[2]; /*y2'=y2*/
return(1);
}
void Euler1(t,y,n,h,k,z)
int n; /*整型变量，微分方程组中方程的个数，也是未知函数的个数*/
int k; /*整型变量。积分步数(包括起始点这一步)*/
double t; /*双精度实型变量,对微分方程进行积分的起始点t0*/
double h; /*双精度实型变量。积分步长*/
double y[]; /*双精度实型一维数组，长度为n。存放n个未知函数yi在起始点t0处的函数值*/
double z[]; /*双精度实型二维数组，体积为nxk。返回k个积分点(包括起始点)上的未知函数值*/
{
extern int Func();
int i,j;
double *d;
d=malloc(n*sizeof(double));
if(d == NULL)
{
printf("内存分配失败\n");
exit(1);
}
/*将方程组的初值赋给数组z[i*k]*/
for (i=0; i<=n-1; i++)
z[i*k]=y[i];
for (j=1; j<=k-1; j++)
{
Func(y,d); /*求出f(x)*/
for(i=0; i<=n-1; i++)
y[i]=z[i*k+j-1]+h*d[i];
Func(y,d);
for (i=0; i<=n-1; i++)
d[i]=z[i*k+j-1]+h*d[i];
for (i=0; i<=n-1; i++)
{
y[i]=(y[i]+d[i])/2.0;
z[i*k+j]=y[i];
}
}
free(d);
return;
}
void Euler2(t,h,y,n,eps)
int n;
double t,h,eps,y[];
{
int i,j,m;
double hh,p,q,*a,*b,*c,*d;
a=malloc(n*sizeof(double));
b=malloc(n*sizeof(double));
c=malloc(n*sizeof(double));
d=malloc(n*sizeof(double));
hh=h;
m=1;
p=1.0+eps;
for (i=0; i<=n-1; i++) a[i]=y[i];
while (p>=eps)
{
for (i=0; i<=n-1; i++)
{
b[i]=y[i];
y[i]=a[i];
}
for (j=0; j<=m-1; j++)
{
for (i=0; i<=n-1; i++)
c[i]=y[i];
Func(y,d);
for (i=0; i<=n-1; i++)
y[i]=c[i]+hh*d[i];
Func(y,d);
for (i=0; i<=n-1; i++)
d[i]=c[i]+hh*d[i];
for (i=0; i<=n-1; i++)
y[i]=(y[i]+d[i])/2.0;
}
p=0.0;
for (i=0; i<=n-1; i++)
{
q=fabs(y[i]-b[i]);
if (q>p) p=q;
}
hh=hh/2.0; m=m+m;
}
free(a);
free(b);
free(c);
free(d);
return;
}
main()
{
int i,j;
double y[3],z[3][11],t,h,x,eps;
y[0]=-1.0; /*初值y0(0)=-1.0*/
y[1]=0.0; /*初值y1(0)=-1.0*/
y[2]=1.0; /*初值y2(0)=-1.0*/
t=0.0; /*起始点t=0*/
h=0.01; /*步长为0.01*/
eps = 0.00001;
Euler1(t,y,3,h,11,z);
printf("定步长欧拉法结果：\n");
for (i=0; i<=10; i++)
{
x=i*h;
printf("t=%5.2f\t ",x);
for(j=0; j<=2; j++)
printf("y(%d)=%e ",j,z[j][i]);
printf("\n");
}
y[0]=-1.0; /*重新赋初值*/
y[1]=0.0;
y[2]=1.0;
printf("变步长欧拉法结果：\n");
printf("t=%5.2f\t ",t);
for (i=0; i<=2; i++)
printf("y(%d)=%e ",i,y[i]);
printf("\n");
for (j=1; j<=10; j++)
{
Euler2(t,h,y,3,eps);
t=t+h;
printf("t=%5.2f\t ",t);
for (i=0; i<=2; i++)
printf("y(%d)=%e ",i,y[i]);
printf("\n");
}
}
trust100@ubuntu:~/test/clanguage$ ./a.out 定步长欧拉法结果： t= 0.00 y(0)=-1.000000e+00 y(1)=0.000000e+00 y(2)=1.000000e+00 t= 0.01 y(0)=-9.999500e-01 y(1)=1.000000e-02 y(2)=9.900500e-01 t= 0.02 y(0)=-9.998000e-01 y(1)=1.999900e-02 y(2)=9.801990e-01 t= 0.03 y(0)=-9.995500e-01 y(1)=2.999600e-02 y(2)=9.704460e-01 t= 0.04 y(0)=-9.992001e-01 y(1)=3.999000e-02 y(2)=9.607901e-01 t= 0.05 y(0)=-9.987502e-01 y(1)=4.998000e-02 y(2)=9.512302e-01 t= 0.06 y(0)=-9.982005e-01 y(1)=5.996501e-02 y(2)=9.417655e-01 t= 0.07 y(0)=-9.975509e-01 y(1)=6.994401e-02 y(2)=9.323949e-01 t= 0.08 y(0)=-9.968016e-01 y(1)=7.991602e-02 y(2)=9.231176e-01 t= 0.09 y(0)=-9.959526e-01 y(1)=8.988004e-02 y(2)=9.139326e-01 t= 0.10 y(0)=-9.950040e-01 y(1)=9.983508e-02 y(2)=9.048389e-01 变步长欧拉法结果： t= 0.00 y(0)=-1.000000e+00 y(1)=0.000000e+00 y(2)=1.000000e+00 t= 0.01 y(0)=-9.999500e-01 y(1)=9.999875e-03 y(2)=9.900499e-01 t= 0.02 y(0)=-9.998000e-01 y(1)=1.999875e-02 y(2)=9.801988e-01 t= 0.03 y(0)=-9.995500e-01 y(1)=2.999563e-02 y(2)=9.704457e-01 t= 0.04 y(0)=-9.992001e-01 y(1)=3.998950e-02 y(2)=9.607896e-01 t= 0.05 y(0)=-9.987503e-01 y(1)=4.997938e-02 y(2)=9.512296e-01 t= 0.06 y(0)=-9.982005e-01 y(1)=5.996426e-02 y(2)=9.417648e-01 t= 0.07 y(0)=-9.975510e-01 y(1)=6.994314e-02 y(2)=9.323941e-01 t= 0.08 y(0)=-9.968017e-01 y(1)=7.991503e-02 y(2)=9.231167e-01 t= 0.09 y(0)=-9.959527e-01 y(1)=8.987892e-02 y(2)=9.139315e-01 t= 0.10 y(0)=-9.950041e-01 y(1)=9.983383e-02 y(2)=9.048378e-01
标签：02,00,01,int,double,步长,实例,90,欧拉 From： https://www.cnblogs.com/mapstar/p/16721715.html