恳请指点:四阶龙格库塔解常微分方程不断改变步长,运行有时候出结果,有时候程序没有结果

d更改很多次,查了不少资料,仍没有结果,恳请大神能帮忙指点一下,或者该参考什么方向的书籍?非常感

谢,当C[0]=200;时

//功能:方程组1和方程组2解不断逼近,t不断改变方程组1,A不断改变方程组2迭代初值,误差满足要求,

输出此时A、t。现在在步长为0.01,C[0]=100;时运行出t=1.1,A=43.0,改变C[0]=200就没有结果了1、参看

资料变了很多步长,最后面的是看的相关论文改的根据误差选取不同的步长,不出结果2,Db[i]和Dc[i]同时

满足,A[i]为负值(理论上为正)
#include
#include
#include
#define PI 3.141593
#define Vx y[1]
#define Vy y[2]
#define X y[3]
#define Y y[4]
#define f1(y1,y2,y3,y4,t) (-0.1109*y1)
#define f2(y1,y2,y3,y4,t) (-0.1109*y2-9.81)//方程2
#define f3(y1,y2,y3,y4,t) (y1)
#define f4(y1,y2,y3,y4,t) (y2)
void main(void)
{
int i,j,V,m=4,n=151;
double t,A[150],Db[150],Dc[150],B[150],C[150],G[150],M[150],r[150],R[150],d,y[150],y0

[150],k1[150],k2[150],k3[150],k4[150];
B[0]=0;
C[0]=100;
G[0]=100;
M[0]=0;
M[1]=0;
M[2]=-10;
V=120;
d=0.01;

t=0;

B[1]=B[0]+M[0]*t;//方程
C[1]=C[0]+M[1]*t;
G[1]=G[0]+M[2]*t;
r[1]=sqrt(B[1]*B[1]+G[1]*G[1]);
R[1]=sqrt(r[1]*r[1]+C[1]*C[1]);
A[1]=atan(C[1]/r[1]);
A[1]=A[1]*180/PI;//弧度化角度
A[1]=cos(A[1]*PI/180);//输入必须为弧度
y0[1]=V*A[1];//初始条件
A[1]=atan(C[1]/r[1]);
A[1]=A[1]*180/PI;//弧度化为度
A[1]=sin(A[1]*PI/180);//输入必须为弧度
y0[2]=V*A[1];
y0[3]=0;
y0[4]=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[i];
Dc[i]=Y-C[i];
//传递t不断更新预测初值
B[i+1]=B[0]+M[0]*t;
C[i+1]=C[0]+M[1]*t;
G[i+1]=G[0]+M[2]*t;
r[i+1]=sqrt(B[i+1]*B[i+1]+G[i+1]*G[i+1]);
R[i+1]=sqrt(r[i+1]*r[i+1]+C[i+1]*C[i+1]);
if((0<Db[i])&&(Db[i]<=1))//误差选择1m只能保证一个误差,另一个相差很大
{
//if((0<Dc[i])&&(Dc[i]<=1))//加上此判断条件时角度为负值,按理论来说不应该是负值
A[i]=atan(Y/X);//此时为弧度
A[i]=A[i]*180/PI;
if(i%1==0)
printf("t=%.1f,角度A=%.1f,方程2X=%.1f,Y=%.1f,方程1r=%.1f,C=%.1f\n",t,A[i],X,Y,r[i

+1],C[i+1]);

}
else//重新更改方程2的迭代初值,
{
A[i+1]=asin(A[i]);
A[i+1]=A[i+1]+Dc[i]/R[i];

    A[i+1]=A[i+1]*180/PI;
    A[i+1]=cos(A[i+1]*PI/180);
    y0[1]=V*A[i+1];
    A[i+1]=asin(A[i]);
    A[i+1]=A[i+1]+Dc[i]/R[i];
    A[i+1]=A[i+1]*180/PI;
    A[i+1]=sin(A[i+1]*PI/180);
    y0[2]=V*A[i+1];
    y0[3]=0;
    y0[4]=0;
    for(j=1;j<=m;j++);
        y[j]=y0[j];
        //  加上下面程序(改变步长d)则运行停止工作,没有时运行良好?
/*  if((0.3*r[i]<X)&&(X<0.6*r[i]))
    {

        d=0.7*d;
    if((0.6*r[i]<X)&&(X<0.9*r[i]))
        d=0.4*d;
    if((0.9*r[i]<X)&&(X<r[i]))
        d=0.1*d;

    }*/




}

}
}

2个回答

还是不显示,include后面是stdio.h、math.h、stdlib.h

不知道为什么一句句粘过来,或者自己重新打也是有些语句发出来就没有了,换了IE浏览器,网页总是错误,试了三四次都掉字母,这次比较少点,前三句是
#include
#include
#include,,恳请大神指点

Csdn user default icon
上传中...
上传图片
插入图片
抄袭、复制答案,以达到刷声望分或其他目的的行为,在CSDN问答是严格禁止的,一经发现立刻封号。是时候展现真正的技术了!
立即提问
相关内容推荐