|
|
发表于 2017-8-1 21:35:10
|
显示全部楼层
确实可以用“中文”实现你所述的数学问题。
比如，我在推算太阳行星运动时，就是用中文加拼音命名的，我觉得效果也不错，还提简捷的，这非线性微分方程组的数值解，这个方程组比较大，有近百个变量。
我用动力学方法计算n体问是，表示方法：
积分表示为 jifen
差分表示为 cafen
行星：XX
外推：waitui
....
我不觉得用中文有多么不好。因为使用C++，所以关键词是英文的，其它都差不多是中文的了。
//---------------------------------------------------------------------------
#pragma hdrstop
#pragma argsused
#include<stdio.h>
#include<conio.h>
#include<math.h>
/**************
1、程序名称：寿星行星积分器
2、版本 V1.1
3、算法来自于《天体力学引论》,作者：牛顿的信徒们
第154页开始,详见：
·纯数值方法
·数值积分的形式
·数值积分的起步问题
·数值积分的累进、性质、应用等章节
4、本程序使用牛顿向前差分公式,而不是使用牛顿中心差分公式,所以算法与原书有所不同
5、本程序包含相对论改正,不含天体形状修正,全部看作质点
6、本程序在C++Builder6中调试通过
7、设计：许剑伟
2009.6.25早晨最后修改
**************/
/**************
====质量与时间的单位====
本程序没有定义质量的单位，单位由初始值决定,但简化了万有与力与加速度关系的
表达,a = G*M/r^2换为a = m/r^2,即G*M换为m, r单位AU,m采用特制单位.当太阳的m取值
0.295912208285591095E-03,则时间单位是"天",九大行星G*M的精密值取自DE405星历的
头文件.
算式a = m/r^2被应用于本程序.表示质量为m的天体作用于相距r的天体产生加速度a.
**************/
/**************
====步长控制====
时元DT(即步长)的大小不同,积分精度也不同,稳定性测试结果如下:
每圈100步,这样1000圈(千年)后角度差10的-7次方弧度,周期误差10的-9
每圈200步,这样1000圈(千年)后角度差10的-10次方弧度,周期误差10的-12
每圈2000步,这样1000圈(千年)后角度差10的-11次方弧度,周期误差10的-12
多步外推的效果与单步外推的效果差不多,但速度要快得多,CN=11时8步外推效果较好
以上精度不是指与实际星历比较的精度
**************/
const CN=11; //最大插值点个数,阶数为CN-1
const PN=10; //天体个数
long double dot(long double *a,long double *b) //矢量点乘
{ return a[0]*b[0] + a[1]*b[1] + a[2]*b[2]; }
//时元设置
//定义二体总质量为1,距离为1,做圆周运动时周期为基准周期,其值为2PI,对应365.2568983263平太阳日
long double DT=2*M_PI/200;
//数值积分系数的计算
long double XSv[CN],XSs[CN]; //v和s的积分系数
void dxss(){//创建多项式表，在程序执行之初调用，CN改变则应dxss应重新调用
if(XSs[0]) return;
int i,j,n;
for(n=0;n<CN;n++){ //建立积分系数
long double jx,a[CN],v=0,s=0,f=0; //(x-f)(x-f-1)...(x-f-n+1)/n!多项式相乘的积分系数,由0积到-1
if(n==0) { XSv[0]=-1;XSs[0]=.5; continue;}
a[n]=1;a[n-1]=f; //计算多项式系数
for(i=n-2;i>=0;i--)
for(j=i,f--,a=0;j<n;j++)
a[j]+=a[j+1]*f;
for(i=1,jx=1;i<=n;i++) jx*=i; //计算n!
for(i=0;i<=n;i++){
if(i%2==0) v-=a/(i+1), s+=a/(i+1)/(i+2);
else v+=a/(i+1), s-=a/(i+1)/(i+2);
}
XSv[n]=v/jx; XSs[n]=s/jx;
}
}
//由等间距插值序列求非插值点的值
long double fx(long double *y,long double t){
int i,j; long double p[CN],v,c;
for(i=0;i<CN;i++) p=y; //复制到差分表
for(i=1,v=p[0],c=1;i<CN;i++){//i是阶数
for(j=0;j<CN-i;j++) //求各阶差分
p[j]=p[j+1]-p[j];
c*=(t-i+1)/i,v+=c*p[0];
}
return v;
}
//天体类
//----------差分与积分--------------
void cafen(long double *p){//建立差分表
int i,j,L;
for(i=1;i<CN;i++){ //i是列
L=(CN-i)*CN;
for(j=0;j<L;j+=CN)//求各阶差分
p[j+i]=p[j+i+CN-1]-p[j+i-1];
}
for(i=1;i<CN;i++){ //补成正方形，i是行
L=i*CN+CN-1; //起始位置从本行末开始
p[L]=p[L-CN];
for(j=1;j<i;j++)
p[L-j]=p[L-j-CN]+p[L-j-CN+1];
}
}
void jifen(long double *p,long double *s,long double *v,int k1=0){ //利用差分表求积分
//s[CN],v[CN]传入积分实始值,s[0至CN-1],v[0至CN-1]传回
int i,j,L;
for(i=k1;i<CN;i++){ //从第k1行开始积分
if(i) v=v[i-1], s=s[i-1];
else v=v[CN], s=s[CN];
for(j=0,L=i*CN;j<CN;j++)
v-=p[L+j]*XSv[j], s-=p[L+j]*XSs[j];
s+=v;
}
}
void waitui(long double *p,long double *s,long double *v,int step){ //外推差分表step步,可以1步,快步最大使用CN-1步(精度不错)
int i,j,c=CN*step,L=(CN-step)*CN;
s[CN]=s[step-1], v[CN]=v[step-1];
for(i=0;i<L;i++) p=p[i+c]; //下移st行
for(i=L,j=0;i<CN*CN;i++,j++){ //按行补成正方形
if(j==CN-1) {j-=CN; p=p[i-CN]; continue; }
p=p[i-CN]+p[i-CN+1];
}
jifen(p,s,v,CN-step); //只需补充积分step行就够
}
class XX{//行星结构
public:
long double cfb[3][CN*CN]; //差分表,0列表示加速度
long double v[3][CN+1],s[3][CN+1]; //最后一个为初态常数,即6个差分表可直接使用的轨道根数
long double Es[CN]; //势能表
long double ba[3][CN]; //备份加速度表
long double m; //天体质量
long double g[6]; //起始轨道根数,g[0-2]位置,g[3-5]速度
long double z[6]; //坐标输出
int bh;
XX(){bh=0;}
long double aBak(){//备份加速度表，并返回当前加速度与原加速度的相对误差(返回精度值)
int i,k; long double max[3]={0.0},c,d;
for(k=0;k<3;k++)
for(i=0;i<CN;i++){
c = cfb[k][i*CN];
d=fabsl(c-ba[k]);
ba[k]=c;
if(max[k]<d) max[k]=d; //最大差值
}
return (max[0]+max[1]+max[2])/sqrtl( ba[0][0]*ba[0][0] + ba[1][0]*ba[1][0] + ba[2][0]*ba[2][0] );
}
void CF() { for(int i=0;i<3;i++) cafen(cfb); } //三维的差分
void JF() { for(int i=0;i<3;i++) jifen(cfb,s,v); } //三维s,v积分
void WT(n){ for(int i=0;i<3;i++) waitui(cfb,s,v,n); } //外推n步
void QB() { //起步外推
//起步外推：把初值s、v复制给各个插值点的s、v,备份表置零,下轮计算a (牛顿定律s算a)
//普通外推:对已积分的表外推，外推后仍是积分后的表,即a推s,下轮计算a(牛顿定律s算a)
//起步外推前调用前，应先初始化g，即初速度各位置
int i,k;
for(k=0;k<3;k++){ //重算加速度后做差分，相当于差分表高阶置零(初始化差分表)
s[k][CN]=g[k]/DT/DT, v[k][CN]=g[k+3]/DT;
for(i=0;i<CN;i++)
s[k]=s[k][CN], v[k]=v[k][CN], ba[k]=0;
}
}
void aCls(int sv=0){//加速度和能量置零,以累加方式重算其它天体对它的加速度前须置零
int i,k;
for(k=0;k<3;k++)
for(i=0;i<CN*CN;i+=CN) cfb[k]=0;
for(i=0;i<CN;i++) Es=0;
}
void svOut(int n){ //转为实际速度及位置
for(int k=0;k<3;k++){
z[k] =s[k][n]*DT*DT;
z[k+3]=v[k][n]*DT;
}
}
};
class TYX{ //太阳系类
public:
XX G[PN]; //九行星
int gn,st; //实际计算天体个数及步长
TYX(){
dxss(); //创建积分系数表
gn=2;
}
void initXX(long double *c,int n){ //行星参数初始化,c质量及根数表,n参与计算的天体实际个数
int i,j,k;
for(i=0;i<n;i++){
j=i*7;
G.m=c[j];
for(k=0;k<6;k++) G.g[k]=c[j+1+k];
}
gn=n;
}
long double aBak(){ //备份加速度并传回误差量
long double w,max=0;
for(int i=0;i<gn;i++){ w=G.aBak(); if(w>max) max=w; }
return max;
}
void aCls(){ for(int i=0;i<gn;i++) G.aCls(); } //加速度置零
void QB() { for(int i=0;i<gn;i++) G.QB(); } //整体起步外推
void WT(n) { for(int i=0;i<gn;i++) G.WT(n); } //整体普通外推
void CF() { for(int i=0;i<gn;i++) G.CF(); } //整体差分
void JF() { for(int i=0;i<gn;i++) G.JF(); } //整体积分
void a_AB(XX &A,XX &B,int hn){ //万有引力产生的加速度,hn为行号
int k;
long double Ad[3],Bd[3],Au[3],Bu[3],Av[3],Bv[3],Aa[3],Ba[3],Ah[3],Bh[3];
long double ABv,Adv,Bdv,Ada,Bda,Av2,Bv2,r,r3;
for(k=0;k<3;k++)
Av[k] = A.v[k][hn] * DT, //A速度矢量
Bv[k] = B.v[k][hn] * DT, //B速度矢量
Aa[k] = A.cfb[k][hn*CN], //A加速度矢
Ba[k] = B.cfb[k][hn*CN], //B加速度矢
Bd[k] = (B.s[k][hn] - A.s[k][hn]) * DT*DT, Ad[k] = -Bd[k], //B相对位置,A相对位置
Bu[k] = Bv[k]-Av[k], Au[k] = -Bu[k], //B相对速度,A相对速度
Ah[k] = 4*Av[k]-3*Bv[k],
Bh[k] = 4*Bv[k]-3*Av[k];
r = dot(Bd,Bd); r = sqrtl(r); r3=r*r*r; //两点间距离
Av2 = dot(Av,Av); //A速度的平方
Bv2 = dot(Bv,Bv); //B速度的平方
ABv = dot(Av,Bv); //A与B速度点乘
Bdv = dot(Bd, Bv); Bdv/=r; Bdv*=Bdv;
Bda = dot(Bd, Ba);
Adv = dot(Ad, Av); Adv/=r; Adv*=Adv;
Ada = dot(Ad, Aa);
long double B1 =4*A.Es[hn] + B.Es[hn] + Av2 + 2*Bv2 - 4*ABv - 1.5*Bdv + Bda/2;
long double A1 =4*B.Es[hn] + A.Es[hn] + Bv2 + 2*Av2 - 4*ABv - 1.5*Adv + Ada/2;
long double B2 = dot(Bd, Ah);
long double A2 = dot(Ad, Bh);
long double B3 = 3.5*r*r, A3 = B3;
long double xz=1;
if(A.bh==0&&B.bh==3||A.bh==3&&B.bh==0){ //太阳与地月系的相互作用力修正
xz=1+0.6377/385/385/80;
}
long double c=299792.458/149597870.691*86400,c2=c*c; //光速
for(k=0;k<3;k++){ //得到加速度
A.cfb[k][hn*CN] += B.m/r3*( Bd[k] + (B1*Bd[k] + B2*Bu[k] + B3*Ba[k])/c2 )*xz;
B.cfb[k][hn*CN] += A.m/r3*( Ad[k] + (A1*Ad[k] + A2*Au[k] + A3*Aa[k])/c2 )*xz;
}
}
void Es_AB(XX &A, XX &B, int hn){ //势能计算
int i,j,k;
long double d[3],r;
for(k=0;k<3;k++) d[k] = (A.s[k][hn] - B.s[k][hn]) * DT*DT; //相对位置矢量
r = sqrtl(dot(d,d)); //两点间距离
A.Es[hn] -= B.m/r;
B.Es[hn] -= A.m/r;
}
void a_All(){ //所有天体间的引力作用
int i,j,hn;
long double r[CN*CN];
for(hn=0;hn<CN;hn++){
for(i=0;i<gn-1;i++) //引力势能计算
for(j=i+1;j<gn;j++) Es_AB(G,G[j],hn);
for(i=0;i<gn-1;i++) //计算所有天体间的力作用，排列组合起来计算
for(j=i+1;j<gn;j++) a_AB(G,G[j],hn);
}
}
void calc(int n,int st){ //积分推步n次,每次st步
//起步外推的步数是CN步,因下次外推时从st步开始
//所以起步外推得到的有效采集点可认为只有st个,其它数据的下次外推时还可取得
int i,j,k;
for(i=0;i<n;i++){
if(!i) QB(); else WT(st); //i==0起步外推,i==1普通外推
for(j=0;j<25;j++){ //整体加速度计算,直到稳定(达预定的精度)。限制次数,实际最多15次就可以稳定
aCls(); a_All(); CF(); JF(); //加速度计算, 差分, 再积分
if(aBak()<1E-13) break;
}
}
}
}S;
void test(){ //二体运动测试
long double cs0[2*7]={
//本组参数为二体运动,周期为2*PI
4, 1,0,0, 0, 1,0, //天体1: 质量(G*M)，x,y,z, vx,vy,vz
4, -1,0,0, 0,-1,0 //天体2: 质量(G*M)，x,y,z, vx,vy,vz
};
DT=2*M_PI/200; //设置步长
int st=8; //每次推步的步数,如:推步200次,则总步数为200*st,如果每圈200步,则积分了st圈
S.initXX(cs0,2); //初始化行星参数
S.calc(200,st); //推步
/*数据输出:
进行n次推步后,共推进n*st步,但差分表只保留最后CN步,编号为0,1...CN-1
如果需要所有的推步的结果,应在calc()及时保存每步结果
差分表中各步与总步数的计数关系
第 0 个对应 (n-1)*st+1 = n*st - (st - 1)
第st-1个对应 n*st
*/
XX *p=S.G;
p[0].svOut(st-1); //取出第0天体第st-1个作坐标,即n*st步坐标
p[1].svOut(st-1); //取出第1天体第st-1个作坐标,即n*st步坐标
printf("%20.17Lf %20.17Lf %20.17Lf\r\n",p[0].z[0],p[0].z[1],p[0].z[2]); //天体0位置print
printf("%20.17Lf %20.17Lf %20.17Lf\r\n",p[1].z[0],p[1].z[1],p[1].z[2]); //天体1位置print
}
//=============九体运动测试===================
void zxsd(long double *cs,int n){ //求参数表中的质心速度
int i,k;
long double *p,v[6]={0},m=0;
for(i=0;i<n;i++){
p=cs+7*i;
for(k=0;k<6;k++) v[k]+=p[0]*p[k+1];
m+=p[0];
}
for(i=0;i<6;i++) v/=m;
printf("SSB速度%20.17Lf %20.17Lf %20.17Lf\r\n",v[3],v[4],v[5]); //天体0置print
}
void test2(){
long double cs[PN*7]={ //支持10体
0.295912208285591095E-03, 0.450250815623389356E-02, 0.767074700932378838E-03, 0.266056805177027132E-03, -0.351748209645188673E-06, 0.517762539958483020E-05, 0.222910185439166520E-05,
0.491254745145081187E-10, 0.361762714603509283E+00, -0.907819677295860494E-01, -0.857149831817633490E-01, 0.336749391398414016E-02, 0.248945204676488743E-01, 0.129463006886503581E-01,
0.724345248616270270E-09, 0.612751941341837969E+00, -0.348365368494968186E+00, -0.195278288980229919E+00, 0.109520683616991937E-01, 0.156176843652620547E-01, 0.633110555360105742E-02,
0.899701134671249882E-09, 0.120517417295400386E+00, -0.925838478939452814E+00, -0.401540220801811010E+00, 0.168112683039926443E-01, 0.174830931338161642E-02, 0.758202868990537761E-03,
0.954953510577925806E-10, -0.110186074282858815E+00, -0.132759945613255570E+01, -0.605889132614203740E+00, 0.144816530597351029E-01, 0.242463117760296197E-03, -0.281520734248005334E-03,
0.282534590952422643E-06, -0.537970689883591202E+01, -0.830480581460151468E+00, -0.224828700228372841E+00, 0.109201154300885245E-02, -0.651811656579268268E-02, -0.282078316536504749E-02,
0.845971518568065874E-07, 0.789439244197905232E+01, 0.459647780162694541E+01, 0.155869757353026683E+01, -0.321755523930330892E-02, 0.433580985895527619E-02, 0.192864656565384452E-02,
0.129202491678196939E-07, -0.182653983068220320E+02, -0.116194450551811457E+01, -0.250103483937377857E+00, 0.221188417417765429E-03, -0.376247593284657683E-02, -0.165101470306800195E-02,
0.152435890078427628E-07, -0.160550425837682305E+02, -0.239421812161785681E+02, -0.940015672354883058E+01, 0.264277104336691632E-02, -0.149831445535920970E-02, -0.679041903017956112E-03,
0.218869976542596968E-11, -0.304833196039993162E+02, -0.872478355495796887E+00, 0.891156304098993246E+01, 0.322210447723588042E-03, -0.314357030215202046E-02, -0.107794882973930419E-02
};
zxsd(cs,10); //显示SSB速度
printf("计算中...\r\n");
DT=0.5; //设置步长,单位天
int i, n=100, st=8; //每次推步的步数st,如:推步200次,则总步数为200*st,如果每圈200步,则积分了st圈
S.initXX(cs,10); //初始化行星参数
for(i=0;i<PN;i++) S.G.bh=i; //给行星编号,太阳必须为0,地球必须为3
S.calc(n,st); //推步
/*数据输出:
进行n次推步后,共推进n*st步,但差分表只保留最后CN步,编号为0,1...CN-1
如果需要所有的推步的结果,应在calc()及时保存每步结果
差分表中各步与总步数的计数关系
第 0 个对应 (n-1)*st+1 = n*st - (st - 1)
第st-1个对应 n*st
*/
XX *p=S.G;
char *mc[10]={"太阳","水星","金星","地月","火星","木星","土星","天王","海王","冥王"};
printf("历元 %Lf\r\n",2451545-11144.5+DT*st*n);
for(i=0;i<PN;i++){
p.svOut(st-1); //取出第i天体第st-1个作坐标,即n*st步坐标
printf("%s %20.15Lf %20.15Lf %20.15Lf\r\n",mc,p.z[0],p.z[1],p.z[2]); //天体置print
}
}
//---------------------------------------------------------------------------
void main(){
test2();
getch();
}
|
|
IP卡
狗仔卡
显身卡
