求助：循环角度不断改变初值，解常微分方程

Question

四阶龙格库塔解微分方程，
#include
#include
#include
#define PI 3.141593
#define Vx y[1]
#define Vy y[2]
#define X y[3]
#define Y y[4]
//方程组2
#define f1(y1,y2,y3,y4,t) (-0.1109*y1)//Vx(Vr)
#define f2(y1,y2,y3,y4,t) (-0.1109*y2-9.81)//Vy（Vh）
#define f3(y1,y2,y3,y4,t) (y1)//X(r)
#define f4(y1,y2,y3,y4,t) (y2)//Y(h)
void main(void)
{
int k,s,i,j,V,m=4,n=150;
double Db[151],Dc[151],T,t,a[151],A[151],A1[151],B[151],C[151],G[151],M

[151],r,R,d,y[151],y0[151],k1[151],k2[151],k3[151],k4[151];

B[0]=0;
C[0]=100;
G[0]=100;
M[0]=0;
M[1]=0;
M[2]=-10;
    V=120;
d=0.1;         
T=0;

for(s=1;s<=100;s++)
{

B[s]=B[s-1]+M[0]*T;//T时刻方程组1
C[s]=C[s-1]+M[1]*T;
G[s]=G[s-1]+M[2]*T;
r=sqrt(B[s]*B[s]+G[s]*G[s]);
R=sqrt(r*r+C[s]*C[s]);
A[s]=atan(C[s]/r);
A[s]=A[s]*180/PI;
A[s]=cos(A[s]*PI/180);//输入必须为弧度
y0[1]=V*A[s];
A[s]=atan(C[s]/r);
A[s]=A[s]*180/PI;//弧度化为度
A[s]=sin(A[s]*PI/180);//输入必须为弧度
y0[2]=V*A[s];
y0[3]=0;
y0[4]=0;
t=0;
for(j=1;j<=m;j++)
y[j]=y0[j];
for(i=1;i<n;i++)
{
for(j=1;j<=m;j++)
{
if(j==1)
k1[j]=d*f1((y[j]),(y[j+1]),(y[j+2]),(y[j+3]),t);
else if(j==2)
k1[j]=d*f2((y[j-1]),(y[j]),(y[j+1]),(y[j+2]),t);
else if(j==3)
k1[j]=d*f3((y[j-2]),(y[j-1]),(y[j]),(y[j+1]),t);
else if(j==4)
k1[j]=d*f4((y[j-3]),(y[j-2]),(y[j-1]),(y[j]),t);
}
for(j=1;j<=m;j++)
{
if(j==1)
k2[j]=d*f1((y[j]+0.5*k1[j]),(y[j+1]+0.5*k1[j+1]),(y[j+2]+0.5*k1[j+2]),(y[j

+3]+0.5*k1[j+3]),(t+0.5*d));
else if(j==2)
k2[j]=d*f2((y[j-1]+0.5*k1[j-1]),(y[j]+0.5*k1[j]),(y[j+1]+0.5*k1[j+1]),(y[j

+2]+0.5*k1[j+2]),(t+0.5*d));
else if(j==3)
k2[j]=d*f3((y[j-2]+0.5*k1[j-2]),(y[j-1]+0.5*k1[j-1]),(y[j]+0.5*k1[j]),(y[j

+1]+0.5*k1[j+1]),(t+0.5*d));
else if(j==4)
k2[j]=d*f4((y[j-3]+0.5*k1[j-3]),(y[j-2]+0.5*k1[j-2]),(y[j-1]+0.5*k1[j-1]),

(y[j]+0.5*k1[j]),(t+0.5*d));
}
for(j=1;j<=m;j++)
{
if(j==1)
k3[j]=d*f1((y[j]+0.5*k2[j]),(y[j+1]+0.5*k2[j+1]),(y[j+2]+0.5*k2[j+2]),(y[j

+3]+0.5*k2[j+3]),(t+0.5*d));
else if(j==2)
k3[j]=d*f2((y[j-1]+0.5*k2[j-1]),(y[j]+0.5*k2[j]),(y[j+1]+0.5*k2[j+1]),(y[j

+2]+0.5*k2[j+2]),(t+0.5*d));
else if(j==3)
k3[j]=d*f3((y[j-2]+0.5*k2[j-2]),(y[j-1]+0.5*k2[j-1]),(y[j]+0.5*k2[j]),(y[j

+1]+0.5*k2[j+1]),(t+0.5*d));
else if(j==4)
k3[j]=d*f4((y[j-3]+0.5*k2[j-3]),(y[j-2]+0.5*k2[j-2]),(y[j-1]+0.5*k2[j-1]),

(y[j]+0.5*k2[j]),(t+0.5*d));
}
for(j=1;j<=m;j++)
{
if(j==1)
k4[j]=d*f1((y[j]+k3[j]),(y[j+1]+k3[j+1]),(y[j+2]+k3[j+2]),(y[j+3]+k3[j+3]),

(t+d));
else if(j==2)
k4[j]=d*f2((y[j-1]+k3[j-1]),(y[j]+k3[j]),(y[j+1]+k3[j+1]),(y[j+2]+k3[j+2]),

(t+d));
else if(j==3)
k4[j]=d*f3((y[j-2]+k3[j-2]),(y[j-1]+k3[j-1]),(y[j]+k3[j]),(y[j+1]+k3[j+1]),

(t+d));
else if(j==4)
k4[j]=d*f4((y[j-3]+k3[j-3]),(y[j-2]+k3[j-2]),(y[j-1]+k3[j-1]),(y[j]+k3[j]),

(t+d));
}
for(j=1;j<=m;j++)
y[j]=y[j]+((k1[j]+2.0*k2[j]+2.0*k3[j]+k4[j])/6.0);
t=i*d;
Db[i]=X-r;
Dc[i]=Y-C[s];//误差选择5m，方程组2的值X、Y需要不断接近方程组1解值C、G,达到精度输出此

时t、A
if(0<Dc[i]<=5)
{
A[i]=asin(A[i]);//此时为弧度
A[i]=A[i]*180/PI;
if(i%10==0)
printf("时间=%.2lf,射角A=%.1lf,击中点X=%.1lf,Y=%.1lf\n",t,A[i],X,Y);
}
else//修正A，更新迭代初值
{
A[s]=atan(C[s]/r);
A[s]=A[s]+Db[i]/R;
A[s]=A[s]*180/PI;
A[s]=cos(A[s]*PI/180);
y0[1]=V*A[s];
A[s]=atan(C[s]/r);//弹丸射角
A[s]=A[s]+Db[i]/R;
A[s]=A[s]*180/PI;
A[s]=sin(A[s]*PI/180);
y0[2]=V*A[s];
y0[3]=0;
y0[4]=0;
for(j=1;j<=m;j++);
y[j]=y0[j];

}
}
//超过最大迭代次数时，重新解方程组1
T=T+n*d;

}
}

Answer 1

运行结果A是乱码，并且怎么好多t从1到最终的循环，不是应该满足精度时才输出t吗