Solve/Plot Second Order DiffEq w/ Two Inital Conditions

I have to solve and plot the result of a second order differential equation using ode45 in MATLAB. I have used ode45 before and cannot figure out how to put a second inital condition in my statements when using a simple code like this (an example of code I have used before):

t2span = [5 10];
v02 = 0; %one initial condition
u2=0; %constant
[t2,v2] = ode45(@(t2,v2) u2-abs(v2)*v2, t2span, v02); %simple ode45 statement

This is the equation given: 𝑥̈ =6sin(𝑡)−𝑥^5 −0.6𝑥̇ , and I must plot x(t) from t=0 to t=240seconds, with initial conditions x(0)=2, 𝑥̇(0)=3.

This is what I have, it works but when I try to do the second part of the problem it makes me think it is incorrect (the second part makes me split the diffeq into its three terms and solve them and then plot them on a combined graph).

M=@(t,X)[X(2);6*sin(t)-X(1).^5-0.6*X(2)] %creating a function
sol = ode45(M,[0 240],[2 3]); %solving using ode45
fplot(@(x)deval(sol,x,1), [0, 240]) %plotting

Please let me know if this is correct or if there is an easier way to do this because I am not sure if this is correct. I got this code from looking up information on ode45 and other ways of solving diffeqs in MATLAB.

Answer 1

[William Rose]

@Erin Hayes,

Your code is good! Instead of using sol on line 2, and instead of fplot() on line 3, you could do as follows. The form below is easier for me to read and understand. I guess that is because it is what I am used to.

Vector t has the time values of the solution. Array x has two columns: [x1,x2]=[x,dx/dt].

M=@(t,X)[X(2);6*sin(t)-X(1).^5-0.6*X(2)];   %create a function
[t,x] = ode45(M,[0 240],[2 3]);             %solve using ode45
plot(t,x(:,1),'-b');                        %plot results
xlabel('Time'); ylabel('x(t)'); grid on

Same result as with your version, of course.

[William Rose]

@Erin Hayes, you said "the second part makes me split the diffeq into its three terms and solve them and then plot them on a combined graph"

Do you mean you want to plot together? Then do this.

M=@(t,X)[X(2);6*sin(t)-X(1).^5-0.6*X(2)];   %create a function
[t,x] = ode45(M,[0 240],[2 3]);             %solve using ode45
%estimate d2x/dt2 at each point using the values of x(:,2)=dx/dt on either side
d2xdt2=(x(3:end,2)-x(1:end-2,2))./(t(3:end)-t(1:end-2));
plot(t,x(:,1),'-r',t,x(:,2),'-g',t(2:end-1),d2xdt2,'-b');  %plot results
xlabel('Time'); ylabel('x(t)'); legend('x','dx/dt','d2x/dt2'); grid on

Try.

[Erin Hayes]

Thank you for your help! Yeah I guess the second part of this question was not too clear the way that I said it. This is what the question is asking and I especially cannot seem to put these three together. I am not sure how to add these types of variables because they are not scalars.

Numerically find the solutions for the following 3 cases. Use initial conditions x(0)=2/3, 𝑥̇(0)=3/3 for all three cases. Plot x1(t)+x2(t)+x3(t).

𝑥1'' = 6sin(𝑡)

𝑥2'' = −𝑥2^5

x3'' = −0.6*𝑥'3

[William Rose]

@Erin Hayes, You have three separate sets of diff eqs to solve. You save the answers. Then you plot the sum of the answers. This will be tricky because the time base for x1 is probably going to be different from the time base for x2 and for x3. Therefore you should save the three solutions, then use interp1() to interpolate each solution to a common time base. Then you can add them up and plot them versus the comon time base. OR you can pass a vector of many times to ode45(), instead of passing tspan. Then ode45() will interpolate its solution to the time vector you passed. Then you won't have to use interp1().

Case 1:

M=@(t,X)[X(2);6*sin(t)];    %create a function
tvec=0:.5:240;              %vector of times
x0=[2/3,1];                 %initial conditions
[t,x] = ode45(M,tvec,x0);   %solve using ode45
x1=x(:,1);                  %save x(:,1) for plotting later

Case 2:

M=@(t,X)[X(2);-X(1).^5];    %create a function
[t,x] = ode45(M,tvec,x0);   %solve using ode45
x2=x(:,1);                  %save x(:,1) for plotting later

Case 3:

M=@(t,X)[X(2);-0.6*X(2)];   %create a function
[t,x] = ode45(M,tvec,x0);   %solve using ode45
x3=x(:,1);                  %save x(:,1) for plotting later

Plot the sum of the results:

plot(t,x1+x2+x3,'-r');

That is not what I expected for the plot. Interesting!

[Erin Hayes]

Thank you! That is what I got as well but I was not sure if it was correct. Thanks so much!

[William Rose]

@Erin Hayes, You're welcome. Good luckwith your work.