计算方法（1~3章）
2012-04-19 21:39
148 查看
写完这些程序还是很有感触的，以前从来没有考虑过误差之类的东西，认为double 16位的精度无所不能~
又想起ac poj 3155那个纠结的时刻，把精度调高居然还wa了！！去问sphinx，他说他也有好多题目是这样，然后我居然在poj 3155的discuss里说这题的测试数据有误。。。- -！汗，直到前几个星期我才意识到错的还是我，学了数值分析才知道有些时候并不是精度越高结果就越准确，在控制double的精度的时候不要让它小于1e-14，在考虑截断误差的同时不要忘了舍入误差，也许poj 3155就是因为在精度设置过高的时候舍入误差产生了很大的影响才wa的吧~
又想起ac poj 3155那个纠结的时刻，把精度调高居然还wa了！！去问sphinx，他说他也有好多题目是这样，然后我居然在poj 3155的discuss里说这题的测试数据有误。。。- -！汗，直到前几个星期我才意识到错的还是我，学了数值分析才知道有些时候并不是精度越高结果就越准确，在控制double的精度的时候不要让它小于1e-14，在考虑截断误差的同时不要忘了舍入误差，也许poj 3155就是因为在精度设置过高的时候舍入误差产生了很大的影响才wa的吧~
/* 拉格朗日插值
*/
#include<iostream>
using namespace std;
const int maxn = 1001;
double x[maxn];
double y[maxn];
double c[maxn]; // ²åֵϵÊýc
double x0,ans;
int n;
double Lagrange()
{
double temp;
for(int k=0;k<=n;k++) //±éÀúËùÓеĵã
{
temp = 1.0;
for(int u=0;u<=n;u++) if(u != k)
{
temp *= x[k] - x[u];
}
c[k] = y[k] / temp;
}
ans = 0.0;
for(int k=0;k<=n;k++)
{
temp = c[k];
for(int u=0;u<=n;u++) if(u != k)
{
temp *= x0 - x[u];
}
ans += temp;
}
return ans;
}
void init() //³õʼ»¯Ã¿¸öµã
{
printf("ÊäÈëÐèÒªÇó½âµÄx0: -> ");
cin>>x0;
printf("\nÇëÊäÈëÒÑÖªµãµÄ¸öÊý: -> ");
cin>>n;
n--;
printf("\n·Ö±ðÊäÈëÿ¸öµã:\n");
for(int i=0;i<=n;i++)
{
printf("\nµÚ%d¸öµã: -> ",i);
cin>>x[i]>>y[i];
}
return ;
}
int main()
{
init();
printf("\nLagrange²åÖµ½á¹ûΪ : ");
cout<<Lagrange()<<endl;
getchar();
getchar();
return 0;
}
/*---- ²âÊÔÊý¾Ý --------*
115
3
100 10
121 11
144 12
----------------------------*//* 埃特金插值
*/
#include<iostream>
#include<iomanip>
using namespace std;
const int maxn = 1001;
double f[maxn][maxn];
double x[maxn];
double y[maxn];
double x0;
int n;
void Aitken()
{
cout.setf(ios::fixed);
cout<<endl<<endl<<f[0][0]<<endl;
for(int k=1;k<=n;k++)
{
cout<<f[0][k]<<" | ";
for(int u=1;u<=k;u++)
{
f[u][k] = (x0 - x[k])*f[u-1][u-1]/(x[u-1]-x[k]) + (x0-x[u-1])*f[u-1][k]/(x[k]-x[u-1]);
if(u!=k)
{
cout<<setprecision(6)<<f[u][k]<<" | ";
}
else
{
cout<<setprecision(6)<<f[u][k]<<endl;
}
}
}
return ;
}
void init() //³õʼ»¯Ã¿¸öµã
{
printf("ÊäÈëÐèÒªÇó½âµÄx0: -> ");
cin>>x0;
printf("\nÇëÊäÈëÒÑÖªµãµÄ¸öÊý: -> ");
cin>>n;
n--;
printf("\n·Ö±ðÊäÈëÿ¸öµã:\n");
for(int i=0;i<=n;i++)
{
printf("\nµÚ%d¸öµã: -> ",i);
cin>>x[i]>>y[i];
}
for(int i=0;i<=n;i++)
{
f[0][i] = y[i];
}
return ;
}
int main()
{
init();
Aitken();
getchar();
getchar();
return 0;
}
/*----- ²âÊÔÊý¾Ý -------*
0.462
5
0.3 0.29850
0.4 0.39646
0.5 0.49311
0.6 0.58813
0.7 0.68122
-----------------------------*//* 牛顿插值
*/
#include<iostream>
using namespace std;
const int maxn = 1001;
double x[maxn];
double y[maxn];
double c[maxn]; //ÐÞÕýϵÊýc
double x0;
int n;
double Newton()
{
double tmp;
double tmp2;
double tmp3;
c[0] = y[0];
for(int k=1;k<=n;k++)
{
tmp = 0.0;
for(int u=0;u<k;u++)
{
tmp2 = c[u];
for(int i=0;i<u;i++)
{
tmp2 *= x[k]-x[i];
}
tmp += tmp2;
}
tmp3 = 1.0;
for(int i=0;i<k;i++)
{
tmp3 *= x[k]-x[i];
}
c[k] = (y[k]-tmp) / tmp3;
}
double ans=0.0;
double temp;
for(int k=0;k<=n;k++)
{
temp = c[k];
for(int i=0;i<k;i++)
{
temp *= x0-x[i];
}
ans += temp;
}
return ans;
}
void init() //³õʼ»¯Ã¿¸öµã
{
printf("ÊäÈëÐèÒªÇó½âµÄx0: -> ");
cin>>x0;
printf("\nÇëÊäÈëÒÑÖªµãµÄ¸öÊý: -> ");
cin>>n;
n--;
printf("\n·Ö±ðÊäÈëÿ¸öµã:\n");
for(int i=0;i<=n;i++)
{
printf("\nµÚ%d¸öµã: -> ",i);
cin>>x[i]>>y[i];
}
return ;
}
int main()
{
init();
printf("\nÅ£¶Ù²åÖµµÄ½á¹ûΪ : ");
cout<<Newton()<<endl;
getchar();
getchar();
return 0;
}
/*---- ²âÊÔÊý¾Ý --------*
115
3
100 10
121 11
144 12
----------------------------*//*---- 复化辛甫生公式--数值积分 ------
----- ¸´»¯ÐÁ¸¦Éú¹«Ê½ ----*/
#include<iostream>
#include<cmath>
using namespace std;
const int maxn = 111;
const double eps = 1e-12;
double x[maxn];
double h;
double f(double x)
{
if(x<eps)
{
return 1.0;
}
return sin(x) / x;
}
double xfs() // 0 --> 1
{
double h = 1.0/100;
for(int i=1;i<=100;i++)
{
x[i] = x[i-1] + h;
}
double ans = 0.0;
double temp;
for(int i=1;i<=100;i++)
{
temp = h/6.0;
temp *= f(x[i-1]) + 4*f((x[i-1]+x[i])/2.0) + f(x[i]);
ans += temp;
printf("%.2lf | %.2lf | %.6lf \n",x[i-1],x[i],temp);
}
return ans;
}
inline void init()
{
memset(x,0,sizeof(x));
return ;
}
int main()
{
init();
printf("answer is : %.12lf\n\n",xfs());
getchar();
return 0;
}/*---- 复化柯特斯公式--数值积分 ------
----- ¸´»¯¿ÂÌØ˹¹«Ê½ ----*/
#include<iostream>
#include<cmath>
using namespace std;
const int maxn = 111;
const double eps = 1e-12;
double x[maxn];
double h;
double f(double x)
{
if(x<eps)
{
return 1.0;
}
return sin(x) / x;
}
double xfs() // 0 --> 1
{
double h = 1.0/100;
for(int i=1;i<=100;i++)
{
x[i] = x[i-1] + h;
}
double ans = 0.0;
double temp;
for(int i=1;i<=100;i++)
{
temp = h/90.0;
temp *= ( + 7.0 * f(x[i-1])
+ 32.0 * f(x[i-1]+h/4.0)
+ 12.0 * f((x[i-1]+x[i])/2.0)
+ 32.0 * f(x[i]-h/4.0)
+ 7.0 * f(x[i])
);
ans += temp;
printf("%.2lf | %.2lf | %.6lf \n",x[i-1],x[i],temp);
}
return ans;
}
inline void init()
{
memset(x,0,sizeof(x));
return ;
}
int main()
{
init();
printf("answer is : %.12lf\n\n",xfs());
getchar();
return 0;
}/* 复化高斯公式
* ÊýÖµ»ý·Ö
*/
#include<iostream>
#include<cmath>
using namespace std;
const double eps = 1e-8;
const int maxn = 100;
double x[maxn+1];
double f(double x)
{
if(x<eps)
{
return 1.0;
}
return sin(x)/x;
}
double gauss()
{
double h = 1.0/maxn;
double cx1 = ( 1.0 - sqrt(1.0/3.0) ) / 2.0;
double cx2 = ( 1.0 + sqrt(1.0/3.0) ) / 2.0;
double x1,x2,temp;
double ans = 0.0;
for(int i=1;i<=maxn;i++)
{
x1 = x[i-1] + (x[i] - x[i-1])*cx1;
x2 = x[i-1] + (x[i] - x[i-1])*cx2;
temp = f(x1) + f(x2);
temp /= 2.0;
temp *= h;
ans += temp;
}
return ans;
}
void init()
{
double h = 1.0 /maxn;
memset(x,0,sizeof(x));
for(int i=1;i<=maxn;i++)
{
x[i] = x[i-1] + h;
}
return ;
}
int main()
{
init();
printf("answer is : %.12lf\n",gauss());
getchar();
return 0;
}/*------ 龙贝格加速算法 ---------
------ Áú±´¸ñ¼ÓËÙËã·¨ --------*/
#include<iostream>
#include<cmath>
using namespace std;
const double eps = 1e-6;
const int maxn = 22;
const int maxstep = 22;
const int maxc = 1<<22;
double x[maxc];
double t[maxn];
double s[maxn];
double c[maxn];
double r[maxn];
double f(double x)
{
if(x<eps)
{
return 1.0;
}
return sin(x)/x;
}
double gets(int id)
{
s[id] = t[id+1]*4.0/3.0 - t[id]/3.0;
return s[id];
}
double getc(int id)
{
c[id] = s[id+1]*16.0/15.0 - s[id]/15.0;
return c[id];
}
double getr(int id)
{
r[id] = c[id+1]*64.0/63.0 - c[id]/63.0;
return r[id];
}
double lbg()
{
double ans=0.0;
double pre = 0;
double now=1000;
int step = 1;
int steps = 0;
int stepc = 0;
int stepr = 0;
double temp;
double h;
printf("| t | s | c | r |\n");
printf("| %.6lf | \n",t[0]);
while( abs(pre - now) > eps || step < maxstep) //¶þ·Ö²½³¤
{
if(stepc >= 1) pre = now;
h = 1.0/(1<<step);
temp = 0.0;
for(int i=1; i <= (1 << step); i++) if(i & 1)
{
temp += f(i * h);
}
temp *= h;
t[step] = t[step - 1] / 2.0 + temp;
printf("| %.6lf | ",t[step]);
if(step >= 1)
{
gets(step - 1);
printf("%.6lf | ",s[step-1]);
steps = step - 1;
}
if(steps >= 1)
{
getc(steps -1);
stepc = steps -1;
printf("%.6lf | ",c[steps-1]);
}
if(stepc >= 1)
{
now = getr(stepc - 1);
stepr = stepc - 1;
printf("%.6lf | ",r[stepc-1]);
}
step ++;
cout<<endl;
}
return now;
}
void init()
{
memset(t,0,sizeof(t));
memset(s,0,sizeof(s));
memset(c,0,sizeof(c));
memset(r,0,sizeof(r));
t[0] = (f(0) + f(1))/2.0;
return ;
}
int main()
{
init();
printf("\nans is : %.12lf \n",lbg());
getchar();
return 0;
}/*数值微分中点方法
* ±ä²½³¤Öе㷽·¨
*/
#include<iostream>
#include<cmath>
using namespace std;
const double eps = 1e-8;
inline double f(double x)
{
return exp(x);
}
double find()
{
double h = 0.8;
double pre=123456789;
double now=987654321;
double a = 1.0;
double l,r;
while(abs(pre-now) > eps)
{
pre = now;
l = a - h;
r = a + h;
l = f(l);
r = f(r);
now = (r - l) / (2.0 * h);
h /= 2.0;
}
return now;
}
int main()
{
printf("answer is : %.12lf \n",find());
getchar();
return 0;
}/* 两步欧拉格式
*/
#include<iostream>
#include<cmath>
using namespace std;
const double eps = 1e-8;
const int maxn = 10;
double x[maxn+1];
double y[maxn+1];
double f(double x,double y)
{
return y - 2.0*x/y;
}
void Two_step_euler_format()
{
double h = 1.0/maxn;
for(int i=2;i<=maxn;i++)
{
y[i] = y[i-2] + 2.0*h*f(x[i-1],y[i-1]);
}
for(int i=1;i<=maxn;i++)
{
cout<<x[i]<<" | ";
printf("%.6lf\n",y[i]);
}
return ;
}
void init()
{
double h = 1.0/maxn;
x[0] = 0.0;
y[0] = 1.0;
for(int i=1;i<=maxn;i++)
{
x[i] = x[i-1] + h;
}
y[1] = y[0] + h * f(x[0] , y[0]);
return ;
}
int main()
{
init();
Two_step_euler_format();
getchar();
return 0;
}/* 改进欧拉格式
* ³£Î¢·Ö·½³ÌµÄ²î·Ö·½·¨
*/
#include<iostream>
#include<cmath>
using namespace std;
const int maxn = 100;
const double eps = 1e-8;
double x[maxn+1];
double y[maxn+1];
double f(double x,double y)
{
return y - 2.0*x/y;
}
void Improved_euler_scheme()
{
double h = 1.0/maxn;
for(int i=1;i<=maxn;i++)
{
y[i] = y[i-1] + h*f(x[i-1],y[i-1]);
y[i] = y[i-1] + h/2.0*(f(x[i-1],y[i-1]) + f(x[i],y[i]));
}
for(int i=1;i<=maxn;i++)
{
cout<<x[i]<<" | ";
printf("%.6lf\n",y[i]);
}
return ;
}
void init()
{
memset(x,0,sizeof(x));
memset(y,0,sizeof(y));
y[0]=1.0;
double h = 1.0/maxn;
for(int i=1;i<=maxn;i++)
{
x[i] = x[i-1] + h;
}
return ;
}
int main()
{
init();
Improved_euler_scheme();
getchar();
return 0;
}/* 四阶龙格库塔经典格式
* ³£Î¢·Ö·½³Ì
*/
#include<iostream>
#include<cmath>
using namespace std;
const double eps = 1e-8;
const int maxn = 10;
double x[maxn+1];
double y[maxn+1];
double f(double x,double y)
{
return y - 2.0*x/y;
}
void start()
{
double k1,k2,k3,k4;
double mid;
double h = 1.0/maxn;
for(int i=1;i<=maxn;i++)
{
mid = (x[i-1]+x[i])/2.0;
k1 = f(x[i-1],y[i-1]);
k2 = f(mid,y[i-1]+h/2.0*k1);
k3 = f(mid,y[i-1]+h/2.0*k2);
k4 = f(x[i],y[i-1]+h*k3);
y[i] = y[i-1] + h/6.0*(k1 + 2*k2 + 2*k3 + k4);
}
for(int i=1;i<=maxn;i++)
{
cout<<x[i]<<" | ";
printf("%.6lf\n",y[i]);
}
return ;
}
void init()
{
double h = 1.0/maxn;
x[0]=0.0;
y[0]=1.0;
for(int i=1;i<=maxn;i++)
{
x[i] = h + x[i-1];
}
return ;
}
int main()
{
init();
start();
getchar();
return 0;
}/* 变步长的龙格库塔格式
* ΢·Ö·½³Ì²î·Ö·½·¨
*/
#include<iostream>
#include<cmath>
using namespace std;
const double eps = 1e-14;
const int maxn = 1<<20;
double x[maxn+1];
double y[maxn+1];
void init(int step);
double f(double x,double y)
{
return y - 2.0*x/y;
}
void start()
{
double k1,k2,k3,k4;
double mid;
int step = 1;
double h;
double pre=123456789;
double now=987654321;
while(abs(pre-now) > eps && step < maxn)
{
init(step);
h = 1.0/step;
pre = now;
for(int i=1;i<=step;i++)
{
mid = (x[i-1]+x[i])/2.0;
k1 = f(x[i-1],y[i-1]);
k2 = f(mid,y[i-1]+h/2.0*k1);
k3 = f(mid,y[i-1]+h/2.0*k2);
k4 = f(x[i],y[i-1]+h*k3);
y[i] = y[i-1] + h/6.0*(k1 + 2*k2 + 2*k3 + k4);
}
now = y[step];
for(int i=1;i<=step;i++)
{
printf("%.6lf | %.12lf\n",x[i],y[i]);
}
printf("Now it has been divided into %d part\n",step);
getchar();
step<<=1;
}
if(abs(now-pre)<eps)
{
printf("exit successful!\n");
}
else
{
printf("error!\n");
}
return ;
}
void init(int step)
{
double h = 1.0/step;
x[0]=0.0;
y[0]=1.0;
for(int i=1;i<=step;i++)
{
x[i] = h + x[i-1];
}
return ;
}
int main()
{
start();
getchar();
return 0;
}/* 四阶亚当姆斯格式
* ΢·Ö·½³Ì²î·Ö·½·¨
*/
#include<iostream>
using namespace std;
const double eps = 1e-8;
const int maxn = 10;
double x[maxn];
double y[maxn];
double f(double x,double y)
{
return y - 2.0*x/y;
}
void start()
{
double k1,k2,k3,k4;
double mid;
double h = 1.0/maxn;
for(int i=1;i<=3;i++) //ËĽ׾µä¸ñʽ
{
mid = (x[i-1]+x[i])/2.0;
k1 = f(x[i-1],y[i-1]);
k2 = f(mid,y[i-1]+h/2.0*k1);
k3 = f(mid,y[i-1]+h/2.0*k2);
k4 = f(x[i],y[i-1]+h*k3);
y[i] = y[i-1] + h/6.0*(k1 + 2*k2 + 2*k3 + k4);
}
for(int i=4;i<=maxn;i++) //Ñǵ±Ä·Ë¹Ô¤±¨-УÕýϵͳ
{
y[i] = y[i-1] + h/24.0 * (
+ 55.0 * f(x[i-1],y[i-1])
- 59.0 * f(x[i-2],y[i-2])
+ 37.0 * f(x[i-3],y[i-3])
- 9.0 * f(x[i-4],y[i-4])
);
y[i] = y[i-1] + h/24.0 * (
+ 9.0 * f(x[i],y[i])
+ 19.0 * f(x[i-1],y[i-1])
- 5.0 * f(x[i-2],y[i-2])
+ f(x[i-3],y[i-3])
);
}
for(int i=1;i<=maxn;i++)
{
cout<<x[i]<<" | ";
printf("%.6lf\n",y[i]);
}
return ;
}
void init()
{
double h = 1.0/maxn;
x[0]=0.0;
y[0]=1.0;
for(int i=1;i<=maxn;i++)
{
x[i] = x[i-1] + h;
}
return ;
}
int main()
{
init();
start();
getchar();
return 0;
}
内容来自用户分享和网络整理，不保证内容的准确性，如有侵权内容，可联系管理员处理 
相关文章推荐
- 最有效的计算2*16的方法——位运算（左移和右移）
- Oracle各种日期计算方法
- php 计算时间差的方法（原创）包括天、小时、分、秒
- 多孔粘性阻力与内部阻力的计算方法
- 编写一个计算长方体体积的C#WindowsForm，包含一个类，一个方法，三个属性
- C语言简单实现计算字符个数的方法
- 卡特兰数的两种计算方法
- 网站并发量的计算方法
- c/c++在windows下获取时间和计算时间差的几种方法总结
- excel计算加权平均方法
- IP和网段的计算方法
- 笨方法学Python 习题 3: 数字和数学计算
- 差集的几种计算方法
- ASP.NET中彩票项目中的计算复式投注的注数的方法
- [导入]c# 统计方法执行时间，计算缓存读取数据方法的命中率
- 圆轨迹弧度计算方法
- 利用算术（加减乘除）的方法实现从1到9计算等于100
- 图像实时处理常用时间间隔计算方法（微秒级）
- PSD-BPA 暂态稳定仿真批处理计算方法的实现
- C/C++日期时间使用方法，struct tm, time_t, clock_t，计算循环所用时间