Plotting Solution Curve on Direction Field

Hello,

I have the second-order differential equation with initial conditions: y'' + 2y' + y = 0, y(-1) = 0, y'(0) = 0

I need to plot the direction field of the solution to the equation and trace the solution curve corresponding to the initial conditions.

I have created the direction field but I'm not sure how to plot the solution curve over the direction field.

Using the plot() function is giving errors and the ezplot() function doesn't seem to represent what the direction field is showing.

Here is what I have so far

% Finds solution to the DE
syms y(x)
Dy = diff(y);
D2y = diff(y,2);
ode = D2y + 2*Dy + y == 0;
ySol1 = dsolve(ode, y(-1)==0); % Solution to DE applying first initial conditions.
ySol2 = dsolve(ode,Dy(0)==0);  % Solution to DE applying second initial conditions.

% Sets up directional field
[x,y]=meshgrid(-4:0.5:4,-4:0.5:4);
u = y;              % x1' = y'
v = - 2*y - x;      % x2' = y'' = - 2*y' - y
u1 = u./sqrt(u.^2+v.^2);
v1 = v./sqrt(u.^2+v.^2);
quiver(x, y, u1, v1, 0.6)
xlabel('x-axis')
ylabel('y-axis')
axis on
axis([-3.5 3.5 -3.5 3.5]);

% Prints the solution curve corresponding to the initial conditions.
hold on
plot(ySol1)
plot(ySol2)
hold off

Any help is greatly appreciated.

Answer 1

[Torsten]

ySol = dsolve(ode, [y(-1)==0,Dy(0)==0]);
dySol = diff(ySol,x);

instead of

ySol1 = dsolve(ode, y(-1)==0); % Solution to DE applying first initial conditions.
ySol2 = dsolve(ode,Dy(0)==0);  % Solution to DE applying second initial conditions.

and

plot(double(subs(ySol,x,-4:0.5:4)),double(subs(dySol,x,-4:0.5:4)))

instead of

plot(ySol1)
plot(ySol2)

[Jordan Stanley]

Hello @Torsten,

Just to understand what you did.

You solved the DE equation with the conditions and assigned it to ySol then took the derivative of that and assigned it to dySol.

Then plotted both ySol and dySol.

Is that correct?

[Torsten]

Yes. This was the task, wasn't it ?

[Jordan Stanley]

Yes, I just wanted to understand the process.

Also, when I replaced my code with your suggestion I get these errors.

Error using sym/subs>normalize (line 240)

Entries in second argument must be scalar.

Error in sym/subs>mupadsubs (line 166)

[X2,Y2,symX,symY] = normalize(X,Y); %#ok

Error in sym/subs (line 154)

G = mupadsubs(F,X,Y);

Do you see anything in my update code that might be causing these.

% Finds solution to the DE
syms y(x)
Dy = diff(y);
D2y = diff(y,2);
ode = D2y + 2*Dy + y == 0;
ySol = dsolve(ode,[y(-1)==0,Dy(0)==0]);
dySol = diff(ySol,x);

% Sets up directional field
[x,y]=meshgrid(-4:0.5:4,-4:0.5:4);
u = y;              % x1' = y'
v = - 2*y - x;      % x2' = y'' = - 2*y' - y
u1 = u./sqrt(u.^2+v.^2);
v1 = v./sqrt(u.^2+v.^2);
quiver(x, y, u1, v1, 0.6)
xlabel('x-axis')
ylabel('y-axis')
axis on
axis([-3.5 3.5 -3.5 3.5]);

% Prints the solution curve corresponding to the initial conditions.
hold on
plot(double(subs(ySol,x,-4:0.5:4)),double(subs(dySol,x,-4:0.5:4)))
hold off

[Torsten]

Then try

syms y(x)
Dy = diff(y);
D2y = diff(y,2);
ode = D2y + 2*Dy + y == 0;
ySol = dsolve(ode,[y(-1)==0,Dy(0)==0]);
dySol = diff(ySol,x);
ySol = matlabFunction(ySol);
dySol = matlabFunction(dySol);
% Sets up directional field
[x,y]=meshgrid(-4:0.5:4,-4:0.5:4);
u = y;              % x1' = y'
v = - 2*y - x;      % x2' = y'' = - 2*y' - y
u1 = u./sqrt(u.^2+v.^2);
v1 = v./sqrt(u.^2+v.^2);
quiver(x, y, u1, v1, 0.6)
xlabel('x-axis')
ylabel('y-axis')
axis on
axis([-3.5 3.5 -3.5 3.5]);
% Prints the solution curve corresponding to the initial conditions.
hold on
plot(ySol(-4:0.5:4),dySol(-4:0.5:4));  
hold off

[Jordan Stanley]

Hello again,

There are less errors now but the only error I get now is,

Not enough input arguments.

Error in symengine>@(C1,x)exp(-x).*(C1+C1.*x)

I have the same code as above.

[Torsten]

I can't run the code, but I'm surprised that ySol should still have a free constant although 2 conditions are given that the solution should fulfill.

What do you get for ySol from the command

ySol = dsolve(ode,[y(-1)==0,Dy(0)==0]);

I solved your equation with paper and pencil and get

ySol(x) = a*exp(-x)*(1+x)

for arbitrary a.

So you can't plot anything here because the solution is not unique.

[Jordan Stanley]

I get an error message that says... "Array indices must be positive integers or logical values.".

[Jordan Stanley]

Interesting, well thank you for your help.

[Torsten]

What I want to say is: You can't trace a solution curve since the solution is not unique.

What you can do is choose a special solution:

syms y(x)
Dy = diff(y);
D2y = diff(y,2);
ode = D2y + 2*Dy + y == 0;
ySol = dsolve(ode,[y(-1)==0,Dy(0)==0]);
dySol = diff(ySol,x);
ySol = @(x) exp(-x).*(1+x);
dySol = @(x) -exp(-x).*x;
% Sets up directional field
[x,y]=meshgrid(-4:0.5:4,-4:0.5:4);
u = y;              % x1' = y'
v = - 2*y - x;      % x2' = y'' = - 2*y' - y
u1 = u./sqrt(u.^2+v.^2);
v1 = v./sqrt(u.^2+v.^2);
quiver(x, y, u1, v1, 0.6)
xlabel('x-axis')
ylabel('y-axis')
axis on
axis([-3.5 3.5 -3.5 3.5]);
% Prints the solution curve corresponding to the initial conditions.
hold on
plot(ySol(-4:0.5:4),dySol(-4:0.5:4));  
hold off

[Jordan Stanley]

Ah okay, I believe that is what I needed.

Thank you very much.

[Torsten]

And what is ySol you get from MATLAB from the line

ySol = dsolve(ode,[y(-1)==0,Dy(0)==0]);

?

[Jordan Stanley]

I still get the error message, "Array indices must be positive integers or logical values."

[Torsten]

And with

syms y(x)
Dy = diff(y,x);
D2y = diff(y,x,2);
ode = D2y + 2*Dy + y == 0;
ySol(x) = dsolve(ode,[y(-1)==0,Dy(0)==0])

What do you get for ySol ?

[Jordan Stanley]

I get...

ySol(x) =

exp(-x)*(C1 + C1*x)

[Torsten]

Then the following code should work (without prescribing ySol as I did before, but the result should be the same):

syms y(x)
Dy = diff(y,x);
D2y = diff(y,x,2);
ode = D2y + 2*Dy + y == 0;
ySol(x) = dsolve(ode,[y(-1)==0,Dy(0)==0])
ySol = subs(ySol,C1,1);
dySol = diff(ySol,x);
ySol = matlabFunction(ySol);
dySol = matlabFunction(dySol);
% Sets up directional field
[x,y]=meshgrid(-4:0.5:4,-4:0.5:4);
u = y;              % x1' = y'
v = - 2*y - x;      % x2' = y'' = - 2*y' - y
u1 = u./sqrt(u.^2+v.^2);
v1 = v./sqrt(u.^2+v.^2);
quiver(x, y, u1, v1, 0.6)
xlabel('x-axis')
ylabel('y-axis')
axis on
axis([-3.5 3.5 -3.5 3.5]);
% Prints the solution curve corresponding to the initial conditions.
hold on
plot(ySol(-4:0.5:4),dySol(-4:0.5:4));  
hold off

[Jordan Stanley]

Thank you very much for the help.