计算物理算法编程

*/ */
程序运行结果： input a,b=0,1↙ s=3.141593
5.5.2 牛顿—柯特斯求积公式 下面给出三到五阶牛顿—柯特斯求积公式。
实现三阶牛顿—柯特斯求积公式的基本步骤如下：
#include <stdio.h> #define N 5 float func(float x) { float y; y=4.0/(1+x*x); return(y); }
（5）若 i N ，i i 1 ，转（4）；否则，转（6）； （6）输出
y
；
（7）结束。
#include <stdio.h> float func(float x,float y) { return(x-y); } float runge_kutta(float x0,float xn,float y0,int N) { float x,y,y1,y2,h,xh; float d1,d2,d3,d4; int i; x=x0; y=y0; h=(xn-x0)/(float)N; /* 计算步长 */ for(i=1;i<=N;i++) /* 标准四阶龙格-库塔公式 */ { xh=x+h/2; d1=func(x,y); d2=func(xh,y+h*d1/2.0); d3=func(xh,y+h*d2/2.0); d4=func(xh,y+h*d3); y=y+h*(d1+2*d2+2*d3+d4)/6.0; x=x0+i*h; } return(y); }
8.2.2 实现欧拉法的基本步骤
（1ห้องสมุดไป่ตู้输入 x 的区间 x 0 , x n ，区间的等分数 N 以及 y 在 x 0 处 的初始值 y 0 ； （2）对 i 0,1,2, N 1 ，计算：
f i f ( xi , y i ) yi yi h f i
x i x i i h
float simpson(float y[],float h) { float s,s1,s2; int i; s1=y[1]; s2=0.0; for(i=2;i<=N-2;i=i+2) { s1+=y[i+1]; /* 计算奇数项的函数值之和 s2+=y[i]; /* 计算偶数项的函数值之和 } s=y[0]+y[N]+4.0*s1+2.0*s2; return(s*h/3.0); } main() { float a,b,h,s,f[N+1]; scanf("%f,%f",&a,&b); h=(b-a)/( float)N; gedianzhi(f,a,h); s=simpson(f,h); printf("s=%f\n",s); }
d 1 f ( x, y)
d 2 f ( x h , y h* d 1 / 2)
d 3 f ( x h , y h* d 2 / 2)
d 4 f ( x h , y h* d 3)
y y h* (d 1 2 * d 2 2 * d 3 d 4)/6
x x 0 i* h
#include <math.h> double func(double x) { double y; y=x*x+sin(x); return(y); } double legendre_gauss(double a,double b,int m,int n) { double h,hx,y,s,dx,x0; int i,k; static double x[]={-0.9061798459,-0.5384693101,0.0000000000, 0.5384693101,0.9061798459}; static double w[]={0.2369268851,0.4786286705,0.5688888889, 0.4786286705,0.2369268851}; dx=(b-a)/(double)m; hx=dx/2; s=0.0; for(k=0;k<m;k++) { x0=a+(double)k*dx+hx; for(i=0;i<n;i++) { y=x0+x[i]*hx; s=s+w[i]*func(y);} } return(s*hx); }
h ( x n - x0 ) / N
（3）输出 y N ； （4）结束。
#include <stdio.h> float func(float x,float y) { return(y-x); }
float euler(float x0,float xn,float y0,int N) { float x,y,h; int i; x=x0; y=y0; h=(xn-x0)/(float)N; /* 计算步长 */ for(i=1;i<=N;i++) /* 欧拉公式 */ { y=y+h*func(x,y); x=x0+i*h; } return(y); }
#include <stdio.h> #define N 16 /* 等分数 */ float func(float x) { float y; y=4.0/(1+x*x); return(y); } void gedianzhi(float y[],float a,float h) { int i; for(i=0;i<=N;i++) y[i]=func(a+i*h); }
(4) 计算
y 0 y N N 1 s h yi 2 i 1
(5) 输出
s 的值；
(6) 结束。
#define N 16 /* 等分数 */ float func(float x) { float y; y=4.0/(1+x*x); return(y); } void gedianzhi(float y[],float a,float h) { int i; for(i=0;i<=N;i++) y[i]=func(a+i*h); } float trapeze(float y[],float h) { float s; int i; s=(y[0]+y[N])/2.0; for(i=1;i<N;i++) s+=y[i]; return(s*h); }
main() { float a,b,h,s,f[N+1]; clrscr(); printf("input a,b="); scanf("%f,%f",&a,&b); h=(b-a)/(float)N; gedianzhi(f,a,h); s=trapeze(f,h); printf("s=%f\n",s); } 程序运行结果： input a,b=0,1↙ s=3.140942
8.4.2 实现标准四阶龙格—库塔法的基本步骤 （1）输入 x 的区间 x0 , x n ，区间的等分数 N 以及 y 在 x 0 处的初始值 y 0 ； （2） h ( x n - x0 ) /N （3） i 1 ； （4）计算：
xh x h / 2
, x x0 ，y y 0 ；
void gedianzhi(float y[],float a,float b,int n) { float h,s; int i; h=(b-a)/(float)n; for(i=0;i<=n;i++) y[i]=func(a+i*h); }
float nc3(float y[],float a,float b,int n) { float h,s,s0,s1; int i; h=(b-a)/(float)n; s0=s1=0.0; for(i=0;i<=n-3;i=i+3) { s0+=y[i]+y[i+3]; s1+=y[i+1]+y[i+2]; } s=s0+3.0*s1; return(s*h*3.0/8.0); } main() { float a,b,h,s,f[N+1]; float n3,n4,n5; printf("input a,b="); scanf("%f,%f",&a,&b); gedianzhi(f,a,b,3); n3=nc3(f,a,b,3); printf("\n 3-nc 4-nc 5-nc\n"); printf("%f %f %f\n",n3,n4,n5); }
5.3.3 实现辛卜生积分法的基本步骤
(1) 输入区间
[a,b] 的端点的 a, b 值以及分割数 N ；
(i 0 ,1, 2, , N ) ；
(2)将区间 [ a,b]等分成 N个小区间，每一个小区间的长度 h ( b a ) / N ；
(3) 计算每一个等分点的函数值 y i f ( a i h )
5.2.3 实现梯形积分法的基本步骤
(1) 输入区间 [ a,b] 的端点 a, b 值以及分割数 N ； (2) 将区间 [ a,b] 等分成 N 个小区间，每一个小区间的 长度 h ( b a ) / N ； (3) 计算每一个等分点的函数值
yi f ( a i h ) (i 0 ,1, 2, , N )
main() { float x0,xn,y0,e; int N; clrscr(); printf("\ninput n:\n "); /* 输入区间等分数 */ scanf("%d",&N); printf("input x0,xn:\n "); /* 输入x的区间 */ scanf("%f%f",&x0,&xn); printf("input y0:\n "); /* 输入x0处的y的值 */ scanf("%f",&y0); e=euler(x0,xn,y0,N); printf("y(%f)=%6.4f",y0,e); }