The fundamental thing you are doing wrong is that you don't have the proper size state vector. For a 1st order equation, the state vector has one element. For a 2nd order equation, the state vector has two elements. If you had two coupled 2nd order equations, the state vector has 2x2 = 4 elements. Etc.
In your case, you have a single 2nd order equation, so your state vector will have two elements. I will keep your same variable names for this discussion. You are using theta(1), theta(2), etc. as your solution states:
theta(1) is the solution at time x(1)
theta(2) is the solution at time x(2)
etc.
But this is inadequate as you only have a single element state at each time point. You need two elements at each time point (one for theta and one for theta_dot) because of your 2nd order equation. I find it convenient to use a matrix to hold the solution with the columns as the states at the time points (but you could easily use rows if you wanted). In this case the states would be:
theta = zeros(2,N);
theta(1:2,1) is the solution at time x(1) for theta and theta_dot
theta(1:2,2) is the solution at time x(2) for theta and theta_dot
etc.
I.e., theta(1,1) is theta at time x(1) and theta(2,1) is theta_dot at time x(1)
So your initialization equations of:
theta(1) = pi/10;
theta1(1) = 0;
would instead be:
theta(1,1) = pi/10;
theta(2,1) = 0;
And your state update equation of:
theta(n+1) = theta(n) + h*(1/3*m*L^2*theta2 + k1*L^2*(1-cos(theta))*sin(theta) + k2*L^2*sin(theta)*cos(theta) - m*g*L/2*cos(theta) + c*L^2*(theta1)*(cos(theta))^2) + theta(n);
would instead be something like:
theta(:,n+1) = theta(:,n) + h*(stuff);
The "stuff" would be a 2x1 vector. The first element of this vector will be theta_dot at the current time, and the second element will be theta_doubledot at the current time. This is close to what you have already, but don't use plain theta without subscripts in the right hand side expression since theta is the ENTIRE solution matrix! Instead, use an expression based on the current solution state which is theta(:,n). (for you to figure out)
And to plot just the "theta" angle part and not the "theta_dot" part, do this:
Also, the following is not correct:
The 5 value is the end time, not the number of time points. So this should be something like this instead:
tstop = 5;
N = round(tstop/h) + 1;
Give all this a shot and come back when you have more problems with your updated code.
Best Answer