计算机仿真(龙格库塔方法)的软件VB设计与实现

计算机仿真（龙格库塔方法）的软件VB设计与实现Dima(0To10)AsSingle,b(0To10)AsSingle,c(0To10)AsSingle,d(0To10)AsSingleDime(0To10)AsSingle,h(0To10)AsSingle,p(0To10)AsSingle,q(0To10)AsSingleDimn1AsByte,n2AsByte'n1表示的是y的阶数,n2表示的是输入函数的阶数DimiAsIntegerDimf(0To2000)AsSingle'f(m)=y,也就是各个时刻的y值DimX0(0To12)AsSingle,X1(0To12)AsSingleDimdtAsSingle,uAsSingle'dt为采样周期,u为输入Dimq1AsSingle,p1AsSingle,h1AsSingle,e1AsSingleDimb0(0To10)AsSingle,bn(-4To10)AsSingleDimmAsIntegerDimmaxAsSinglePrivateSubCombo1_Click()n1=Combo1.ListIndex+1'n1表示的是y的阶数Fori=n1+1To8Text1(i).Visible=FalseLabel1(i).Visible=FalseLabel9(i).Visible=FalseNextFori=0Ton1Text1(i).Visible=TrueLabel1(i).Visible=TrueLabel9(i).Visible=TrueNextCombo2.ClearCombo2.Text=请选择输出u的最大阶数Fori=0Ton1Combo2.AddItem((i)&阶)NextEndSubPrivateSubCombo2_Click()n2=Combo2.ListIndex'n2表示的是输入函数的阶数Fori=n2+1To8Text2(i).Visible=FalseLabel2(i).Visible=FalseLabel10(i).Visible=FalseNextFori=0Ton2Text2(i).Visible=TrueLabel2(i).Visible=TrueLabel10(i).Visible=TrueNextEndSubPrivateSubCalculate_Click()dt=Val(Text3.Text)'dt是采样周期u=Val(Text4.Text)n1=Combo1.ListIndex+1Ifn1=0Then'对付忘记选择最高阶数时的情况CallMsgBox(请选择阶数！,48,未选择阶数)ExitSubEndIfFori=0To12'设X0的初值都是0，且i大于8，这一点很重要！X0(i)=0NextFori=0Ton1-1a(i)=-Val(Text1(i).Text)/Val(Text1(n1).Text)'a(i)b(i)c(i)d(i)等等都是计算的中间变量NextFori=0Ton2b0(i)=Val(Text2(i).Text)/Val(Text1(n1).Text)NextIfCombo3.ListIndex=1Thenn1=n1+1Fori=n1-1To1Step-1a(i)=a(i-1)Nexta(0)=0EndIfFori=-4To-1bn(i)=0Nextb(0)=a(0)*a(n1-1)Fori=1Ton1-1b(i)=a(i-1)+a(i)*a(n1-1)Nextc(0)=a(0)*b(n1-1)Fori=1Ton1-1c(i)=b(i-1)+a(i)*b(n1-1)Nextd(0)=a(0)*c(n1-1)Fori=1Ton1-1d(i)=c(i-1)+a(i)*c(n1-1)NextFori=0Ton1-1e(i)=dt*a(i)+dt*dt/2*b(i)+dt*dt*dt/6*c(i)+dt*dt*dt*dt/24*d(i)h(i)=dt*dt/2*a(i)+dt*dt*dt/6*b(i)+dt*dt*dt*dt/24*c(i)p(i)=dt*dt*dt/6*a(i)+dt*dt*dt*dt/24*b(i)q(i)=dt*dt*dt*dt/24*a(i)NextIfn2=0Then'n2为0和不为0的公式是不一样的，表示输入函数u为n2阶'b0(0)=Val(Text2(0).Text)/Val(Text1(n1).Text)m=1Doq1=0:p1=0:h1=0:e1=0Fori=0Ton1-1'这些个也是中间变量q1=q1+q(i)*X0(i+1)p1=p1+p(i)*X0(i+1)h1=h1+h(i)*X0(i+1)e1=e1+e(i)*X0(i+1)NextFori=1Ton1X1(i)=X0(i)+dt*X0(i+1)+dt*dt/2*X0(i+2)+dt*dt*dt/6*X0(i+3)+dt*dt*dt*dt/24*X0(i+4)NextX1(n1)=X1(n1)+e1+(dt+dt*dt/2*a(n1-1)+dt*dt*dt/6*b(n1-1)+dt*dt*dt*dt/24*c(n1-1))*b0(0)*uIfn1=4ThenX1(n1-3)=X1(n1-3)+q1+dt*dt*dt*dt/24*b0(0)*uX1(n1-2)=X1(n1-2)+p1+(dt*dt*dt/6+dt*dt*dt*dt/24*a(n1-1))*b0(0)*uX1(n1-1)=X1(n1-1)+h1+(dt*dt/2+dt*dt*dt/6*a(n1-1)+dt*dt*dt*dt/24*b(n1-1))*b0(0)*uElseIfn1=3ThenX1(n1-2)=X1(n1-2)+p1+(dt*dt*dt/6+dt*dt*dt*dt/24*a(n1-1))*b0(0)*uX1(n1-1)=X1(n1-1)+h1+(dt*dt/2+dt*dt*dt/6*a(n1-1)+dt*dt*dt*dt/24*b(n1-1))*b0(0)*uElseIfn1=2ThenX1(n1-1)=X1(n1-1)+h1+(dt*dt/2+dt*dt*dt/6*a(n1-1)+dt*dt*dt*dt/24*b(n1-1))*b0(0)*uEndIff(m)=X1(1)Fori=1Ton1X0(i)=X1(i)Nextm=m+1LoopUntilm=2000'循环计算2000次ElseFori=(n2+1)To10b0(i)=0NextFori=0Ton1-1bn(i)=b0(i)+b0(n1)*a(i)Nextm=1DoFori=1Ton1X1(i)=X0(i)+dt*X0(i+1)+dt*dt/2*X0(i+2)+dt*dt*dt/6*X0(i+3)+dt*dt*dt*dt/24*X0(i+4)+(h(n1-i)*bn(n1-1)+p(n1-i)*bn(n1-2)+q(n1-i)*bn(n1-3))*uNextFori=1Ton1X1(i)=X1(i)+e(n1-i)*X0(1)+h(n1-i)*X0(2)+p(n1-i)*X0(3)+q(n1-i)*X0(4)+(dt*bn(n1-i)+dt*dt/2*bn(n1-1-i)+dt*dt*dt/6*bn(n1-2-i)+dt*dt*dt*dt/24*bn(n1-3-i))*uNextf(m)=X1(1)Fori=1Ton1X0(i)=X1(i)Nextm=m+1LoopUntilm=2000EndIfmax=f(1)Form=2To1999IfAbs(f(m))Abs(max)Then'开始绘图max=f(m)EndIfNext'Ifmax0ThenPicture1.ClsPicture1.Scale(-20*dt,1.2*max)-(2020*dt,-0.05*max)Picture1.Line(0,0)-(2002*dt,0)Picture1.Line(0,0)-(0,1.18*max)Picture1.Line(2002*dt,0)-(1960*dt,-0.025*max)Picture1.Line(2002*dt,0)-(1960*dt,0.025*max)Picture1.Line(0,1.18*max)-(10*dt,1.1*max)Picture1.Line(0,1.18*max)-(-10*dt,1.1*max)'Else'Picture1.Cls'Picture1.Scale(-5*dt,0.2*Abs(max))-(2002*dt,1.2*max)'Picture1.Line(0,0)-(2002*dt,0)'Picture1.Line(0,0)-(0,1.2*max)'Picture1.Line(2002*dt,0)-(1900*dt,0.1*max)'Picture1.Line(0,1.2*max)-(50*dt,1.1*max)'EndIfForm=2To1999Picture1.Line((m-1)*dt,f(m-1))-(m*dt,f(m)),vbRedNext'If(f(1999)-f(1800))0.00001Then'Picture1.Line(0,f(1999))-(2000*dt,f(1999))'Picture1.CurrentX=0'Picture1.CurrentY=f(1999)+0.1*max'Picture1.Printy∞=&f(1999)'EndIfIfCombo3.ListIndex=0ThenPicture1.Line(0,f(1999))-(2000*dt,f(1999))Picture1.CurrentX=0Picture1.CurrentY=f(1999)+0.1*maxPicture1.Printy∞=&f(1999)ElsePicture1.Line(0,0)-(2000*dt,u*2000*dt)EndIfEndSubPrivateSubCombo3_Click()Text4.Visible=TrueLabel12.Visible=TrueIfCombo3.ListIndex=0ThenLabel5.Caption=请输入阶跃信号u的系数：ElseLabel5.Caption=请输入斜坡信号u的系数：EndIfEndSubPrivateSubdtxnzb_Click()Dimy00AsSingleDimthigemaAsSingleDimtsAsSingleDimessAsSingleIfCombo3.ListIndex=0Theny00=f(1999)thigema=max-y00m=2000Dom=m-1LoopUntilAbs(f(m)-y00)=(0.05*y00)ts=m*dtess=u-y00Frame1.Visible=TrueLabel6.Caption=超调量σ:&thigemaLabel7.Caption=调节时间ts:&tsLabel8.Caption=稳态误差ess:&essElseFrame1.Visible=TrueLabel6.Caption=稳态误差ess:&(u*1999*dt-f(1999))/(u*1999*dt)&%EndIfEndSub