zhang3_hm2.zip_matlab画RK方程_simulink 振动_微分方程的数值求解-四阶龙格库塔方法_振动

x1(1)=0; x2(1)=1; t(1)=0; m=1;k=2;h=0.2; for i=1:20 k11(i)=-k/m*x2(i)^3+1/m*sin(2*t(i)); k12(i)=x1(i); k21(i)=-k/m*(x2(i)+h/2*k11(i))^3+1/m*sin(2*(t(i)+h/2)); k22(i)=x1(i)+h/2*k12(i); k31(i)=-k/m*(x2(i)+h/2*k21(i))^3+1/m*sin(2*(t(i)+h/2)); k32(i)=x1(i)+h/2*k22(i); k41(i)=-k/m*(x2(i)+h*k31(i))^3+1/m*sin(2*(t(i)+h)); k42(i)=x1(i)+h*k32(i); x1(i+1)=x1(i)+h/6*(k11(i)+2*k21(i)+2*k31(i)+k41(i)); x2(i+1)=x2(i)+h/6*(k12(i)+2*k22(i)+2*k32(i)+k42(i)); t(i+1)=t(i)+0.2; end plot(t,x1,'g',t,x2,'r'); legend('dx','x'); xlabel('t'); ylabel('x'); xlim([0,4]);