辛算法

计算物理辛算法和保结构算法是我国著名数学家冯康及其学派在80年代中期系统提出、并完善和发展起来的。采用辛算法和C程序语言实现如下的实际问题。222(1)()2()()0aaadkkFUFd（1）22)1()]([2)(U222()()()0aadFrFd（2）其中，222aUVVW，22222213244ddkdWVVVddd1()exp()VD15.0,137D当0,,1acrk时，（cr为1()exp()VD的点）。方程的解写为)(sin)(1)(2/1AFa(3)但是往往需要()aF的值。因0,,1acrk时此要求采用辛算法和C程序语言求解方程（2）给出一组,()aF的值。提示：方程（2）转化为2()()()()()()aaaaadGrFBrFddFGd初始条件为100aaFcrGcr其中c为常数，由于波函数的归一化，这里取1c，22122k采用数值格式为1111111212112122212323223233323113433343(),,(),,(),,(),.nnnnnnnvGchBFFdhvdhvvchBdhvdhvvchBdhvdhGvchBFdhv1/31(22),12,123420,,,ccccc1/2,d2()/2,d3241,dddd.C语言程序：（在C-free上运行的）#includestdio.h#includemath.hdoubleblta=1.0/137,k=1,D=5.0,e=0.2,p=1e-2;doubleB(doublea){doubleB;B=-2*(e-U(a))+k*(k+1)/(a*a);returnB;}doubleU(doublea){doubleU;U=V(a)-blta*blta/2*((V(a)-e)*(V(a)-e)-W(a));returnU;}doubleV(doublea){doubleV;V=-1/a*exp(-a/D);returnV;}doubleW(doublea){doubleW;W=1/(4*E(a))*(dv2_dp(a)+3*blta*blta/4/E(a)*dv_dp(a)*dv_dp(a)-2*k/a*dv_dp(a));returnW;}doubleE(doublea){doubleE;E=1+blta*blta/2*(1-V(a));returnE;}doubledv_dp(doublea){doubledv_dp;dv_dp=1/a/a*exp(-a/D)+1/a/D*exp(-a/D);returndv_dp;}doubledv2_dp(doublea){doubledv2_dp;dv2_dp=(-2/a/a/a-2/a/a/D-1/a/D/D)*exp(-a/D);returndv2_dp;}doubledeidai(doublep,doubleG,doubleF){doubleh,v1,v2,v3,v4,u1,u2,u3,t,t1,t2,t3,a,b;FILE*fp;inti=1;h=1e-2;a=1.0/(2-pow(2,1.0/3));b=1-2*a;fp=fopen(F:\\xinsuanfa.txt,wr);fprintf(fp,p\t\tG\t\tF\n);fprintf(fp,%5.6lf\t%5.6lf\t%5.6lf\n,p,G,F);for(i=1;;i++){v1=G;u1=F+a/2*v1*h;t1=p+a/2*h;v2=v1-a*h*B(t1)*u1;u2=u1+(a+b)/2*h*v2;t2=t1+(a+b)/2*h;v3=v2-b*h*B(t2)*u2;u3=u2+(a+b)/2*h*v3;t3=t2+(a+b)/2*h;G=v3-a*h*B(t3)*u3;F=u3+a/2*h*v3;p=p+h;printf(p=%5.6lf,G=%5.6lf,F=%5.6lf\n,p,G,F);if(!fp){printf(filecannotbeopened);exit(1);}fprintf(fp,%5.6lf\t%5.6lf\t%5.6lf\n,p,G,F);if(p1.0){break;}}}main(){doublet,G,F;t=(1+2*sqrt(k*k-blta*blta))/2;G=t*pow(p,(t-1));F=pow(p,t);printf(p=%5.6lf,G=%5.6lf,F=%5.6lf\n,p,G,F);deidai(p,G,F);}