Create .m file

---------------------------------------------

function?theta=danbai(t,X)

x=X(1);

dx=X(2);

ddx=-sin(x);

theta=[dx;ddx];

----------------------------------------------

Command window input

>>? [t,Y]=ode45(@danbai,[0?6],[pi/3?-1/2]);

>>?plot(t,Y(:,1),'ro-',t,Y(:,2),'bv-');

>>?legend('\theta-t','d\ theta-t')?

Self-coding Runge Kutta

------------------------------------------------

function?[y,z]=Runge_kutta(a,b,y0,z0,h)

x=a:h:b;< /p>

y(1)=y0;

z(1)=z0;

n=(b-a)/h+1;

for?i=2:n

K(1,1)=f1(x(i-1),y(i-1),z(i-1));

K(2,1)=f2(x(i-1),y(i -1),z(i-1));

K(1,2)=f1(x(i-1)+h/2,y(i-1)+K(1,1)*h/2,z(i-1)+K(2,1)*h/2);

K(2,2)=f2(x(i-1)+h/2,y(i-1)+K(1,1)*h/2,z(i-1)+ K(2,1)*h/2);

K(1,3)=f1(x(i-1)+h/2,y(i-1)+K(1,2)*h/2,z(i-1)+K(2,2)*h/2);

K(2,3)=f2(x(i-1)+h/2,y(i-1)+K(1,2)*h/2,z(i-1)+ K(2,2)*h/2);

K(1,4)=f1(x(i-1)+h,y(i-1)+K(1,3)*h,z(i-1)+K(2,3)*h);

K(2,4)=f2(x(i-1)+h,y(i-1)+K(1,3)*h,z(i-1)+K(2,3)*h).

y(i)=y(i-1)+h/6*(K(1,1)+2*K(1,2)+2*K(1,3)+K(1,4));

z(i)=z(i-1)+h/6*(K(2,1)+2*K(2,2)+2*K(2,3)+K(2,4));

end

y(2)

plot(y,'r')? %θ-t plot

hold?on

plot(f1(x,y,z),'g')? %?dθ-t plot

hold?on

plot(f2(x,y,z),'b-')%d2θ-t plot

%f1.m

function?f1=f1(x,y,z)

f1=z;

%f2.m

function?f2=f2(x,y,z)

f2=-sin(y);

-----------------------------------------

Runge_kutta(0,6,pi/3,-1/2,0.02)