- 技术人生
- 定积分
辛普生求积分：
#include
<
stdio
.
h
>
#include
<math.
h
>
double
Function
(
double
);
double
SIMP1
(
double
,
double
,
int
);
double
SIMP2
(
double
,
double
,
double
);
void
main
()
{
double
a1
,
b1
,
eps
;
int
n1
;
double
Result1
;
double
Result2
;
a1
= 0.0;
b1
= 0.8;
n1
= 4;
eps
= 5e-7;
Result1
=
SIMP1
(
a1
,
b1
,
n1
);
Result2
=
SIMP2
(
a1
,
b1
,
eps
);
printf
(
"
利用定步长辛普生积分结果为：
/n"
);
printf
(
"I1 = %.10f/n"
,
Result1
);
printf
(
"
利用变步长辛普生积分结果为：
/n"
);
printf
(
"I2 = %e/n"
,
Result2
);
}
double
SIMP1
(
a
,
b
,
n
)
double
a
;
double
b
;
int
n
;
{
int
i
;
double
h
,
s
;
h
=(
a
-
b
)/(2*
n
);
s
=0.5*(
Function
(
a
)-
Function
(
b
));
for
(
i
=1;
i
<=
n
;
i
++)
s
+=2*
Function
(
a
+(2*
i
-1)*
h
)+
Function
(
a
+2*
i
*
h
);
return
((
b
-
a
)*
s
/(3*
n
));
}
double
SIMP2
(
a
,
b
,
eps
)
double
a
;
double
b
;
double
eps
;
{
int
k
,
n
;
double
h
,
t1
,
t2
,
s1
,
s2
,
p
,
x
;
n
=1;
h
=
b
-
a
;
t1
=
h
*(
Function
(
a
)+
Function
(
b
))/2;
s1
=
t1
;
while
(1)
{
p
= 0;
for
(
k
=0;
k
<=
n
;
k
++)
{
x
=
a
+(
k
+0.5)*
h
;
p
+=
Function
(
x
);
}
t2
=(
t1
+
h
*
p
)/2;
s2
=(4*
t2
-
t1
)/3;
if
(
fabs
(
s2
-
s1
)>=
eps
)
{
t1
=
t2
;
n
=
n
+
n
;
h
=
h
/2;
s1
=
s2
;
continue
;
}
break
;
}
return
(
s2
);
}
double
Function
(
double
x
)
{
return
(
cos
(
x
));
}
欧拉方程法：
#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"
);
}
}
