?? integral.cs
字號(hào):
/*
* 計(jì)算數(shù)值積分的類(lèi) Integral
*
* 周長(zhǎng)發(fā)編制
*/
using System;
namespace CSharpAlgorithm.Algorithm
{
/**
* 計(jì)算數(shù)值積分的類(lèi) Integral
*
* @author 周長(zhǎng)發(fā)
* @version 1.0
*/
public abstract class Integral
{
/**
* 抽象函數(shù)：計(jì)算積分函數(shù)值，必須在派生類(lèi)中覆蓋該函數(shù)
*
* @param x - 函數(shù)變量
* @return double型，對(duì)應(yīng)的函數(shù)值
*/
public abstract double Func(double x);
/**
* 基本構(gòu)造函數(shù)
*/
public Integral()
{
}
/**
* 變步長(zhǎng)梯形求積法
*
* 調(diào)用時(shí)，須覆蓋計(jì)算函數(shù)f(x)值的虛函數(shù)double Func(double x)
*
* @param a - 積分下限
* @param b - 積分上限，要求b>a
* @param eps - 積分精度要求
* @return double 型，積分值
*/
public double GetValueTrapezia(double a, double b, double eps)
{
int n,k;
double fa,fb,h,t1,p,s,x,t=0;
// 積分區(qū)間端點(diǎn)的函數(shù)值
fa=Func(a);
fb=Func(b);
// 迭代初值
n=1;
h=b-a;
t1=h*(fa+fb)/2.0;
p=eps+1.0;
// 迭代計(jì)算
while (p>=eps)
{
s=0.0;
for (k=0;k<=n-1;k++)
{
x=a+(k+0.5)*h;
s=s+Func(x);
}
t=(t1+h*s)/2.0;
p=Math.Abs(t1-t);
t1=t;
n=n+n;
h=h/2.0;
}
return(t);
}
/**
* 變步長(zhǎng)辛卜生求積法
*
* 調(diào)用時(shí)，須覆蓋計(jì)算函數(shù)f(x)值的虛函數(shù)double Func(double x)
*
* @param a - 積分下限
* @param b - 積分上限，要求b>a
* @param eps - 積分精度要求
* @return double 型，積分值
*/
public double GetValueSimpson(double a, double b, double eps)
{
int n,k;
double h,t1,t2,s1,s2=0,ep,p,x;
// 計(jì)算初值
n=1;
h=b-a;
t1=h*(Func(a)+Func(b))/2.0;
s1=t1;
ep=eps+1.0;
// 迭代計(jì)算
while (ep>=eps)
{
p=0.0;
for (k=0;k<=n-1;k++)
{
x=a+(k+0.5)*h;
p=p+Func(x);
}
t2=(t1+h*p)/2.0;
s2=(4.0*t2-t1)/3.0;
ep=Math.Abs(s2-s1);
t1=t2; s1=s2; n=n+n; h=h/2.0;
}
return(s2);
}
/**
* 自適應(yīng)梯形求積法
*
* 調(diào)用時(shí)，須覆蓋計(jì)算函數(shù)f(x)值的虛函數(shù)double Func(double x)
*
* @param a - 積分下限
* @param b - 積分上限，要求b>a
* @param d - 對(duì)積分區(qū)間進(jìn)行分割的最小步長(zhǎng)，當(dāng)子區(qū)間的寬度
* 小于d時(shí)，即使沒(méi)有滿足精度要求，也不再往下進(jìn)行分割
* @param eps - 積分精度要求
* @return double 型，積分值
*/
public double GetValueATrapezia(double a, double b, double d, double eps)
{
double h,f0,f1,t0,z;
double[] t = new double[2];
// 迭代初值
h=b-a;
t[0]=0.0;
f0=Func(a);
f1=Func(b);
t0=h*(f0+f1)/2.0;
// 遞歸計(jì)算
ppp(a,b,h,f0,f1,t0,eps,d,t);
z=t[0];
return(z);
}
/**
* 內(nèi)部函數(shù)
*/
private void ppp(double x0, double x1, double h, double f0, double f1, double t0, double eps, double d, double[] t)
{
double x,f,t1,t2,p,g,eps1;
x=x0+h/2.0;
f=Func(x);
t1=h*(f0+f)/4.0;
t2=h*(f+f1)/4.0;
p=Math.Abs(t0-(t1+t2));
if ((p<eps)||(h/2.0<d))
{
t[0]=t[0]+(t1+t2);
return;
}
else
{
g=h/2.0;
eps1=eps/1.4;
// 遞歸
ppp(x0,x,g,f0,f,t1,eps1,d,t);
ppp(x,x1,g,f,f1,t2,eps1,d,t);
return;
}
}
/**
* 龍貝格求積法
*
* 調(diào)用時(shí)，須覆蓋計(jì)算函數(shù)f(x)值的虛函數(shù)double Func(double x)
*
* @param a - 積分下限
* @param b - 積分上限，要求b>a
* @param eps - 積分精度要求
* @return double 型，積分值
*/
public double GetValueRomberg(double a, double b, double eps)
{
int m,n,i,k;
double h,ep,p,x,s,q=0;
double[] y = new double[10];
// 迭代初值
h=b-a;
y[0]=h*(Func(a)+Func(b))/2.0;
m=1;
n=1;
ep=eps+1.0;
// 迭代計(jì)算
while ((ep>=eps)&&(m<=9))
{
p=0.0;
for (i=0;i<=n-1;i++)
{
x=a+(i+0.5)*h;
p=p+Func(x);
}
p=(y[0]+h*p)/2.0;
s=1.0;
for (k=1;k<=m;k++)
{
s=4.0*s;
q=(s*p-y[k-1])/(s-1.0);
y[k-1]=p; p=q;
}
ep=Math.Abs(q-y[m-1]);
m=m+1;
y[m-1]=q;
n=n+n;
h=h/2.0;
}
return(q);
}
/**
* 計(jì)算一維積分的連分式法
*
* 調(diào)用時(shí)，須覆蓋計(jì)算函數(shù)f(x)值的虛函數(shù)double Func(double x)
*
* @param a - 積分下限
* @param b - 積分上限，要求b>a
* @param eps - 積分精度要求
* @return double 型，積分值
*/
public double GetValuePq(double a, double b, double eps)
{
int m,n,k,l,j;
double hh,t1,s1,ep,s,x,t2,g=0;
double[] h = new double[10];
double[] bb = new double[10];
// 迭代初值
m=1;
n=1;
hh=b-a;
h[0]=hh;
t1=hh*(Func(a)+Func(b))/2.0;
s1=t1;
bb[0]=s1;
ep=1.0+eps;
// 迭代計(jì)算
while ((ep>=eps)&&(m<=9))
{
s=0.0;
for (k=0;k<=n-1;k++)
{
x=a+(k+0.5)*hh;
s=s+Func(x);
}
t2=(t1+hh*s)/2.0;
m=m+1;
h[m-1]=h[m-2]/2.0;
g=t2;
l=0;
j=2;
while ((l==0)&&(j<=m))
{
s=g-bb[j-2];
if (Math.Abs(s)+1.0==1.0)
l=1;
else
g=(h[m-1]-h[j-2])/s;
j=j+1;
}
bb[m-1]=g;
if (l!=0)
bb[m-1]=1.0e+35;
g=bb[m-1];
for (j=m;j>=2;j--)
g=bb[j-2]-h[j-2]/g;
ep=Math.Abs(g-s1);
s1=g;
t1=t2;
hh=hh/2.0;
n=n+n;
}
return(g);
}
/**
* 高振蕩函數(shù)求積法
*
* 調(diào)用時(shí)，須覆蓋計(jì)算函數(shù)f(x)值的虛函數(shù)double Func(double x)
*
* @param a - 積分下限
* @param b - 積分上限，要求b>a
* @param m - 被積函數(shù)中振蕩函數(shù)的角頻率
* @param n - 給定積分區(qū)間兩端點(diǎn)上的導(dǎo)數(shù)最高階數(shù)＋1
* @param fa - 一維數(shù)組，長(zhǎng)度為n，存放f(x)在積分區(qū)間端點(diǎn)x=a處的各階導(dǎo)數(shù)值
* @param fb - 一維數(shù)組，長(zhǎng)度為n，存放f(x)在積分區(qū)間端點(diǎn)x=b處的各階導(dǎo)數(shù)值
* @param s - 一維數(shù)組，長(zhǎng)度為2，其中s(1)返回f(x)cos(mx)在積分區(qū)間的積分值，
* s(2) 返回f(x)sin(mx)在積分區(qū)間的積分值
* @return double 型，積分值
*/
public double GetValuePart(double a, double b, int m, int n, double[] fa, double[] fb, double[] s)
{
int mm,k,j;
double sma,smb,cma,cmb;
double[] sa = new double[4];
double[] sb = new double[4];
double[] ca = new double[4];
double[] cb = new double[4];
// 三角函數(shù)值
sma=Math.Sin(m*a);
smb=Math.Sin(m*b);
cma=Math.Cos(m*a);
cmb=Math.Cos(m*b);
// 迭代初值
sa[0]=sma;
sa[1]=cma;
sa[2]=-sma;
sa[3]=-cma;
sb[0]=smb;
sb[1]=cmb;
sb[2]=-smb;
sb[3]=-cmb;
ca[0]=cma;
ca[1]=-sma;
ca[2]=-cma;
ca[3]=sma;
cb[0]=cmb;
cb[1]=-smb;
cb[2]=-cmb;
cb[3]=smb;
s[0]=0.0;
s[1]=0.0;
mm=1;
// 循環(huán)迭代
for (k=0;k<=n-1;k++)
{
j=k;
while (j>=4)
j=j-4;
mm=mm*m;
s[0]=s[0]+(fb[k]*sb[j]-fa[k]*sa[j])/(1.0*mm);
s[1]=s[1]+(fb[k]*cb[j]-fa[k]*ca[j])/(1.0*mm);
}
s[1]=-s[1];
return s[0];
}
/**
* 勒讓德－高斯求積法
*
* 調(diào)用時(shí)，須覆蓋計(jì)算函數(shù)f(x)值的虛函數(shù)double Func(double x)
*
* @param a - 積分下限
* @param b - 積分上限，要求b>a
* @param eps - 積分精度要求
* @return double 型，積分值
*/
public double GetValueLegdGauss(double a, double b, double eps)
{
int m,i,j;
double s,p,ep,h,aa,bb,w,x,g=0;
// 勒讓德－高斯求積系數(shù)
double[] t={-0.9061798459,-0.5384693101,0.0,
0.5384693101,0.9061798459};
double[] c={0.2369268851,0.4786286705,0.5688888889,
0.4786286705,0.2369268851};
// 迭代初值
m=1;
h=b-a;
s=Math.Abs(0.001*h);
p=1.0e+35;
ep=eps+1.0;
// 迭代計(jì)算
while ((ep>=eps)&&(Math.Abs(h)>s))
{
g=0.0;
for (i=1;i<=m;i++)
{
aa=a+(i-1.0)*h;
bb=a+i*h;
w=0.0;
for (j=0;j<=4;j++)
{
x=((bb-aa)*t[j]+(bb+aa))/2.0;
w=w+Func(x)*c[j];
}
g=g+w;
}
g=g*h/2.0;
ep=Math.Abs(g-p)/(1.0+Math.Abs(g));
p=g;
m=m+1;
h=(b-a)/m;
}
return(g);
}
/**
* 拉蓋爾－高斯求積法
*
* 調(diào)用時(shí)，須覆蓋計(jì)算函數(shù)f(x)值的虛函數(shù)double Func(double x)
*
* @return double 型，積分值
*/
public double GetValueLgreGauss()
{
int i;
double x,g;
// 拉蓋爾－高斯求積系數(shù)
double[] t={0.26355990, 1.41340290, 3.59642600, 7.08580990, 12.64080000};
double[] c={0.6790941054, 1.638487956, 2.769426772, 4.315944000, 7.104896230};
// 循環(huán)計(jì)算
g=0.0;
for (i=0; i<=4; i++)
{
x=t[i];
g=g+c[i]*Func(x);
}
return(g);
}
/**
* 埃爾米特－高斯求積法
*
* 調(diào)用時(shí)，須覆蓋計(jì)算函數(shù)f(x)值的虛函數(shù)double Func(double x)
*
* @return double 型，積分值
*/
public double GetValueHermiteGauss()
{
int i;
double x,g;
// 埃爾米特－高斯求積系數(shù)
double[] t={-2.02018200, -0.95857190, 0.0,0.95857190, 2.02018200};
double[] c={1.181469599, 0.9865791417, 0.9453089237, 0.9865791417, 1.181469599};
// 循環(huán)計(jì)算
g=0.0;
for (i=0; i<=4; i++)
{
x=t[i];
g=g+c[i]*Func(x);
}
return(g);
}
}
}?? 快捷鍵說(shuō)明
復(fù)制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號(hào)
Ctrl + =
減小字號(hào)
Ctrl + -