岑蓉络阳 发表于 2017-8-1 20:40
从头看到尾，没有发现计算方法专业的朋友。打个比方，计算一个N个未知量的线性方程组或者一个常系数微分方 ...
确实可以用“中文”实现你所述的数学问题。
比如，我在推算太阳行星运动时，就是用中文加拼音命名的，我觉得效果也不错，还提简捷的，这非线性微分方程组的数值解，这个方程组比较大，有近百个变量。
我用动力学方法计算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*b + a*b + a*b; }
//时元设置
//定义二体总质量为1,距离为1,做圆周运动时周期为基准周期,其值为2PI,对应365.2568983263平太阳日
long double DT=2*M_PI/200;
//数值积分系数的计算
long double XSv,XSs; //v和s的积分系数
void dxss(){//创建多项式表，在程序执行之初调用，CN改变则应dxss应重新调用
if(XSs) return;
int i,j,n;
for(n=0;n<CN;n++){ //建立积分系数
long double jx,a,v=0,s=0,f=0; //(x-f)(x-f-1)...(x-f-n+1)/n!多项式相乘的积分系数,由0积到-1
if(n==0) { XSv=-1;XSs=.5; continue;}
a=1;a=f;//计算多项式系数
for(i=n-2;i>=0;i--)
for(j=i,f--,a=0;j<n;j++)
a+=a*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=v/jx;XSs=s/jx;
}
}
//由等间距插值序列求非插值点的值
long double fx(long double *y,long double t){
int i,j; long double p,v,c;
for(i=0;i<CN;i++) p=y; //复制到差分表
for(i=1,v=p,c=1;i<CN;i++){//i是阶数
for(j=0;j<CN-i;j++) //求各阶差分
p=p-p;
c*=(t-i+1)/i,v+=c*p;
}
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=p-p;
}
for(i=1;i<CN;i++){ //补成正方形，i是行
L=i*CN+CN-1; //起始位置从本行末开始
p=p;
for(j=1;j<i;j++)
p=p+p;
}
}
void jifen(long double *p,long double *s,long double *v,int k1=0){ //利用差分表求积分
//s,v传入积分实始值,s,v传回
int i,j,L;
for(i=k1;i<CN;i++){ //从第k1行开始积分
if(i) v=v, s=s;
elsev=v,s=s;
for(j=0,L=i*CN;j<CN;j++)
v-=p*XSv, s-=p*XSs;
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=s, v=v;
for(i=0;i<L;i++) p=p; //下移st行
for(i=L,j=0;i<CN*CN;i++,j++){ //按行补成正方形
if(j==CN-1) {j-=CN; p=p; continue; }
p=p+p;
}
jifen(p,s,v,CN-step); //只需补充积分step行就够
}
class XX{//行星结构
public:
long double cfb; //差分表,0列表示加速度
long double v,s; //最后一个为初态常数,即6个差分表可直接使用的轨道根数
long double Es; //势能表
long double ba; //备份加速度表
long double m; //天体质量
long double g; //起始轨道根数,g位置,g速度
long double z; //坐标输出
int bh;
XX(){bh=0;}
long double aBak(){//备份加速度表，并返回当前加速度与原加速度的相对误差(返回精度值)
int i,k; long double max={0.0},c,d;
for(k=0;k<3;k++)
for(i=0;i<CN;i++){
c = cfb;
d=fabsl(c-ba);
ba=c;
if(max<d) max=d; //最大差值
}
return (max+max+max)/sqrtl( ba*ba + ba*ba + ba*ba );
}
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=g/DT/DT, v=g/DT;
for(i=0;i<CN;i++)
s=s, v=v, ba=0;
}
}
void aCls(int sv=0){//加速度和能量置零,以累加方式重算其它天体对它的加速度前须置零
int i,k;
for(k=0;k<3;k++)
for(i=0;i<CN*CN;i+=CN) cfb=0;
for(i=0;i<CN;i++) Es=0;
}
void svOut(int n){ //转为实际速度及位置
for(int k=0;k<3;k++){
z=s*DT*DT;
z=v*DT;
}
}
};
class TYX{ //太阳系类
public:
XX G; //九行星
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;
for(k=0;k<6;k++) G.g=c;
}
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,Bd,Au,Bu,Av,Bv,Aa,Ba,Ah,Bh;
long double ABv,Adv,Bdv,Ada,Bda,Av2,Bv2,r,r3;
for(k=0;k<3;k++)
Av = A.v * DT, //A速度矢量
Bv = B.v * DT, //B速度矢量
Aa = A.cfb, //A加速度矢
Ba = B.cfb, //B加速度矢
Bd = (B.s - A.s) * DT*DT, Ad = -Bd, //B相对位置,A相对位置
Bu = Bv-Av, Au = -Bu, //B相对速度,A相对速度
Ah = 4*Av-3*Bv,
Bh = 4*Bv-3*Av;
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 + B.Es + Av2 + 2*Bv2 - 4*ABv - 1.5*Bdv + Bda/2;
long double A1 =4*B.Es + A.Es + 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 += B.m/r3*( Bd + (B1*Bd + B2*Bu + B3*Ba)/c2 )*xz;
B.cfb += A.m/r3*( Ad + (A1*Ad + A2*Au + A3*Aa)/c2 )*xz;
}
}
void Es_AB(XX &A, XX &B, int hn){ //势能计算
int i,j,k;
long double d,r;
for(k=0;k<3;k++) d = (A.s - B.s) * DT*DT; //相对位置矢量
r = sqrtl(dot(d,d)); //两点间距离
A.Es -= B.m/r;
B.Es -= A.m/r;
}
void a_All(){ //所有天体间的引力作用
int i,j,hn;
long double r;
for(hn=0;hn<CN;hn++){
for(i=0;i<gn-1;i++) //引力势能计算
for(j=i+1;j<gn;j++) Es_AB(G,G,hn);
for(i=0;i<gn-1;i++) //计算所有天体间的力作用，排列组合起来计算
for(j=i+1;j<gn;j++) a_AB(G,G,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*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.svOut(st-1); //取出第0天体第st-1个作坐标,即n*st步坐标
p.svOut(st-1); //取出第1天体第st-1个作坐标,即n*st步坐标
printf("%20.17Lf %20.17Lf %20.17Lf\r\n",p.z,p.z,p.z); //天体0位置print
printf("%20.17Lf %20.17Lf %20.17Lf\r\n",p.z,p.z,p.z); //天体1位置print
}
//=============九体运动测试===================
void zxsd(long double *cs,int n){ //求参数表中的质心速度
int i,k;
long double *p,v={0},m=0;
for(i=0;i<n;i++){
p=cs+7*i;
for(k=0;k<6;k++) v+=p*p;
m+=p;
}
for(i=0;i<6;i++) v/=m;
printf("SSB速度%20.17Lf %20.17Lf %20.17Lf\r\n",v,v,v); //天体0置print
}
void test2(){
long double cs={ //支持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={"太阳","水星","金星","地月","火星","木星","土星","天王","海王","冥王"};
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,p.z,p.z); //天体置print
}
}
//---------------------------------------------------------------------------
void main(){
test2();
getch();
}
这是力学公式，是微分方程：爱因斯坦－因非尔德－霍夫曼公式。
因为方程复杂，所以我把它写成上述形式，以方便编程。
也就是说，数学公式，也不一定适合编程，最好重新表达为适合计算机的，不然容易错。
用了拼音，最好配上中文注释，这样就好更解了。
xjw01 发表于 2017-8-1 21:43
用了拼音，最好配上中文注释，这样就好更解了。
如果变量、函数、数据结构定义等等直接用中文，很多注释就可以省掉，程序本身更具有可读性，中国人写程序就象写小说一样了，势必大大提高中国人编写复杂软件的能力。
上述的 爱因斯坦－因非尔德－霍夫曼公式，写在纸上，足足写1页纸，即有积分，又有微分。
最后，我还是用中文方式重新表达了这个公式。
比如，和Ea表示势能，配上文字说明，写程序时，就把所有天体势能取和后代入就可以了，而不是表达为一推西格马公式
ssckb 发表于 2017-7-25 20:02
在我记忆中，应该有十多年了，中国人就有人开发了中文版的编程语言，记得好像叫易语言。
不过也没有如楼主 ...
到不是不兼容，VC，vb的都可以调用问题是，编程的难点根本不在那几个字母上，就是你换成了汉字不会还是不会。易语言用的也很多，功能也很强大，很多外挂，病毒都易语言开发的。
yngz 发表于 2017-8-1 22:01
如果变量、函数、数据结构定义等等直接用中文，很多注释就可以省掉，程序本身更具有可读性，中国人写程序 ...
不不不，那是你想当然了，程序代替不了注释，注释不是对程序的描述，而是当时你解决问题的思路过程。
2549608436 发表于 2017-7-25 22:30
你拿汉语注释一下大跃进,人民公社和文化大革命试试?那是已经完成了的程序.
由最高程序员设计,一整套的团队 ...
编程语言就是编程语言，和自然语言没有关系，既没有汉语编程也没有英语编程。
2549608436 发表于 2017-7-25 22:30
你拿汉语注释一下大跃进,人民公社和文化大革命试试?那是已经完成了的程序.
由最高程序员设计,一整套的团队 ...
不懂，不要放屁。我们这是技术论坛，你有你们的地方。
有可以用汉语的开发环境的，单片机和pc的都有，网上搜一下。
多了解，多学习。
比高级语言还在高级的“超高级语言”
好像有中文编程的哦，易语言好像