Solving two second order ODEs

Hi there!

I am trying to solve for U(t) and V(t) for the two second order ODEs

and ,

where a and w are constants and with the initial conditions and .

Then I want to plot the solutions against for a given time interval.

syms U(t) V(t)
%Constants definition
a = 1;
w = 100;
dU=diff(U,t);
dV=diff(V,t);
%Initial Conditions
y0 = [0 0 1 1];
eq1 = diff(U, t, 2) == -a*w*dV;
eq2 = diff(V,t, 2) == a*w*dU;
vars = [U(t); V(t)]
OTV = odeToVectorField([eq1,eq2])
M = matlabFunction(OTV,'vars', {'t','Y'});
interval = [0 5];  %time interval    
ySol = ode45(M,interval,y0);
tValues = linspace(interval(1),interval(2),1000);
yValues1 = deval(ySol,tValues,1); %U(t) solution
yValues2 = deval(ySol,tValues,2); %V(t) solution
plot(yValues1,yValues2)  

Does the above code correctly solve the system of differential equations and initial conditions and plot V(t) against U(t)?

If not, then please let me know what is incorrect. Also plese let me know if there is another way to solve the system of ODEs.

Thank you for your help. Much appreciated.

Answer 1

[Star Strider]

Essentially, yes.

However it could be made a bit more efficient —

syms U(t) V(t)
%Constants definition
a = 1;
w = 100;
dU=diff(U,t);
dV=diff(V,t);
%Initial Conditions
y0 = [0 0 1 1];
eq1 = diff(U, t, 2) == -a*w*dV;
eq2 = diff(V,t, 2) == a*w*dU;
vars = [U(t); V(t)]

vars = 

[OTV,Subs] = odeToVectorField([eq1,eq2])

OTV = 

Subs = 

M = matlabFunction(OTV,'vars', {'t','Y'});
interval = [0 5];  %time interval    
ySol = ode45(M,interval,y0);
tValues = linspace(interval(1),interval(2),1000);
yValues1 = deval(ySol,tValues,1); %U(t) solution
yValues2 = deval(ySol,tValues,2); %V(t) solution
plot(yValues1,yValues2)  

% ---------- Slightly More Efficient: Solves Directly & Avoids The 'deval' Calls ----------
interval = [0 5];  %time interval    
tValues = linspace(interval(1),interval(2),1000);
[t,y] = ode45(M,tValues,y0);

yValues1 = y(:,1); %U(t) solution
yValues2 = y(:,2); %V(t) solution

figure
plot(yValues1,yValues2)
xlabel('$V(t)$', 'Interpreter','latex')
ylabel('$\frac{dV(t)}{dt}$', 'Interpreter','latex')

.

[Jake Barlow]

Hi Star Strider, thank you very much for the answer. Are the axes labels correct in the last plot? Should it say U(t) on the x-axis and V(t) on the y-axis instead?

[Star Strider]

My pleasure!

It was not clear what the first plot represented.

That is where the ‘Subs’ output of odeToVectorField helps, and the reason I always include it in the requested outputs from odeToVEctorVield. According to it, the first two columns (in the second integration that I added) are V and . The U and results would be ‘y(:3)’ and ‘y(:,4)’ in order.

To plot V as a function of U would be —

syms U(t) V(t)
%Constants definition
a = 1;
w = 100;
dU=diff(U,t);
dV=diff(V,t);
%Initial Conditions
y0 = [0 0 1 1];
eq1 = diff(U, t, 2) == -a*w*dV;
eq2 = diff(V,t, 2) == a*w*dU;

[OTV,Subs] = odeToVectorField([eq1,eq2])

OTV = 

Subs = 

M = matlabFunction(OTV,'vars', {'t','Y'});

% ---------- Slightly More Efficient: Solves Directly & Avoids The 'deval' Calls ----------
interval = [0 5];  %time interval    
tValues = linspace(interval(1),interval(2),1000);
[t,y] = ode45(M,tValues,y0);
yValues1 = y(:,3); %U(t) solution
yValues2 = y(:,1); %V(t) solution
figure
plot(yValues1,yValues2)
xlabel('$U(t)$', 'Interpreter','latex')
ylabel('$V(t)$', 'Interpreter','latex')

Check that to be certain, however it appears to be the desired result to me.

.

[Jake Barlow]

Hi Star Strider, thank you very much for your comment. It is clearer now!

[Star Strider]

My pleasure!

Answer 2

[Paul]

I think there are a few mistakes in the code

syms U(t) V(t)
%Constants definition
a = 1;
w = 100;
dU=diff(U,t);
dV=diff(V,t);
%Initial Conditions
y0 = [0 0 1 1];
eq1 = diff(U, t, 2) == -a*w*dV;
eq2 = diff(V,t, 2) == a*w*dU;
vars = [U(t); V(t)]

vars = 

[OTV,S] = odeToVectorField([eq1,eq2]);
S

S = 

Note that S is ordered [V dV U dU], so that should be the ordering of the solution of ode45

M = matlabFunction(OTV,'vars', {'t','Y'});
interval = [0 5];  %time interval
% note the IC's are in the same order as S
ySol = ode45(M,interval,[0 1 0 1],odeset('MaxStep',0.0001,'InitialStep',0.0001));
tValues = linspace(interval(1),interval(2),10000);
yValues1 = deval(ySol,tValues,1); % V(t) solution
yValues2 = deval(ySol,tValues,2); % Vdot(t) solution
yValues3 = deval(ySol,tValues,3); % U(t) solution
yValues4 = deval(ySol,tValues,4); % Udot(t) solution
figure
plot(yValues3,yValues1)  % plot U vs V
xlabel('U');ylabel('V')

An exact solution can be computed using dsolve()

sol = dsolve([eq1; eq2],[U(0)==0; V(0)==0; dU(0)==1; dV(0)==1])

sol = struct with fields:
    V: sin(100*t)/100 - cos(100*t)/100 + 1/100
    U: cos(100*t)/100 + sin(100*t)/100 - 1/100

Ufunc = matlabFunction(sol.U);
Vfunc = matlabFunction(sol.V);
plot(Ufunc(tValues),Vfunc(tValues))
xlabel('U');ylabel('V')

% compare numerical and exact solutions
figure
plot(tValues,yValues3-Ufunc(tValues),tValues,yValues1-Vfunc(tValues))

[Jake Barlow]

Hi Paul, thank you very much for your answer.