Example 1. (System of Linear Equations) To find the general solution to the system of equations
x
0 = 2x + 3y
y
0 = 2x + y,
type at the command line
syms t x y
[x y]
= dsolve(’Dx = 2x+3y’, ’Dy = 2x+y’)
and Matlab returns the following:
x =
(3C2exp(4t))/2 – C1exp(-t)
y =
C1exp(-t) + C2exp(4t)
Note that this is the general solution:
X~ =
x
y
= c1
−1
1
e
−t + c2
3/2
1
e
4t
.
To solve the same equation with the initial condition x(0) = 4, y(0) = 1, type the following two lines and the
output follows:
[x y]
= dsolve(’Dx = 2x+3y’, ’Dy = 2x+y’, ’x(0)=4’, ’y(0)=1’)
x =
exp(-t) + 3exp(4t)
y =
2exp(4*t) – exp(-t)
And this is the particular solution:
X~ =
x
y
−1
1
e
−t +
3
2
e
4t
.
The solutions of x, y (0 ≤ t ≤ 10) and its phase portrait are plotted:
syms x(t) y(t)
[x y]
= dsolve(’Dx = 2x+3y’, ’Dy = 2x+y’, ’x(0)=4’, ’y(0)=1’);
T=[0:0.01:10]’;
X = double( subs(x, t, T) );
Y = double( subs(y, t, T) );
figure(1)
plot(T,X,’r-’,T,Y,’k-’,’linewidth’,2)
legend(’x(t)’,’y(t)’)
xlabel(’t’,’fontsize’,16)
ylabel(’x and y’,’fontsize’,16)
MATH 430 Lab 2 Fall 2019
figure(2)
plot(X,Y,’r-’,’linewidth’,2)
title(’Phase Portrait’, ’fontsize’,16);
xlabel(’x’,’fontsize’,16)
ylabel(’y’,’fontsize’,16)
Figure 1: Plot of Example 1. Left are solutions. Right is the phase portrait.
Example 2. (Numerical Solutions of a linear system) To find the general solution to the system of equations
x
0 = 2x + 3y
y
0 = 2x + y,
in range 0 ≤ t ≤ 5, with the initial condition x(0) = 4, y(0) = 1, with increment h = 0.1. To solve this system
numerically, type at the command lines
f1 = @(t,x,y) 2x + 3y
f1 =
function_handle with value:
@(t,x,y)2x+3y
f2 = @(t,x,y) 2x + y
f2 =
function_handle with value:
@(t,x,y)2x+y
t0 = 0; x0 = 4; y0 = 1; t1 = 5;
h = 0.1; n = (t1-t0)/h;
t = zeros(n+1,1);
x = zeros(n+1,1);
y = zeros(n+1,1);
t(1) = t0; x(1) = x0; y(1) = y0;
for i=1:n
x(i+1) = x(i) + hf1(t(i),x(i),y(i));
y(i+1) = y(i) + h*f2(t(i),x(i),y(i));
t(i+1) = t(i) + h;
end
figure(3)
plot(t,x,’r-’,t,y,’k-’,’linewidth’,2)
2
MATH 430 Lab 2 Fall 2019
legend(’x(t)’,’y(t)’)
xlabel(’t’,’fontsize’,16)
ylabel(’x and y’,’fontsize’,16)
The solutions are plotted in the figure below. We can also plot its phase portrait:
figure(4)
plot(x,y,’r-’,’linewidth’,2)
title(’Phase Portrait’, ’fontsize’,16);
xlabel(’x’,’fontsize’,16)
ylabel(’y’,’fontsize’,16)
Figure 2: Plot of Example 2. Left are solutions. Right is the phase portrait.
Example 3. (Eigenpairs of a linear system) For a linear system
x
0 = 2x + 3y
y
0 = 2x + y.
Matlab can calculate the eigenvalues and eigenvectors:
A = [2 3; 2 1]
A =
2 3
2 1
[v lambda]
= eig(A)
v =
0.8321 -0.7071
0.5547 0.7071
lambda =
4 0
0 -1
Note that it returns two eigenpairs:
λ1 = 4, v1 =
0.8321
0.5547
, λ2 = −1, v2 =
−0.7071
0.7071
.
Example 4. (Matrix exponential) For a given matrix, Matlab can calculate its matrix exponential:
3
MATH 430 Lab 2 Fall 2019
A = [2 3; 2 1]
A =
2 3
2 1
expm(A)
ans =
32.9060 32.5382
21.6921 22.0600
exp(A)
ans =
7.3891 20.0855
7.3891 2.7183
syms t;
B = [1 t; 0 1]
B =
[ 1, t]
[ 0, 1]
expm(B)
ans =
[ exp(1), texp(1)]
[ 0, exp(1)]
exp(B)
ans =
[ exp(1), exp(t)]
[ 1, exp(1)]
Note that in Matlab, “exp” is different from “expm”. “exp” is element-wise operation:
exp a b
c d returns
e
a
e
b
e
c
e
d
while, “expm” is the matrix exponential defined in class:
expm a b
c d returns e
a b
c d
Example 5. (Variation of Constant Formula) For the following linear system
x
0 =
2t 0
0 3
x +
t
1
, x(0) =
0
2
,
4
MATH 430 Lab 2 Fall 2019
First, use the previous method. Type the following in Matlab:
syms t x y
[x y]
= dsolve(’Dx=2tx+t’, ’Dy=3y+1’, ’x(0)=0’, ’y(0)=2’)
x =
exp(t^2)/2 – 1/2
y =
-exp(3t)(exp(-3t)/3 – 7/3)
simplify([x;y]) % show a simplified result
ans =
exp(t^2)/2 – 1/2
(7exp(3*t))/3 – 1/3
Note this is the solution:
x
y
e
t
2
/2 − 1/2
7e
3t/3 − 1/3
Now use the variation of constant formula. Type in Matlab:
syms t;
A = [2t 0; 0 3];
b = [t;1];
t0 = 0;
x0 = [0;2];
Phi = expm(int(A,t,t0,t)); % Here use e^{integral_t0^t A(s) ds} as fundamental matrix
inv_Phi = inv(Phi); % inverse of fundamental matrix
inv_Phi_0 = subs(inv_Phi,t,t0); % value of this inverse when t=t0
x = Phiinv_Phi_0x0 + Phiint(inv_Phib, t, t0, t); % variation of constant formula
x = simplify(x) % show a simplified result
x =
exp(t^2)/2 – 1/2
(7exp(3*t))/3 – 1/3
Lab 2 Assignment:
- Try the above examples. Include the scripts, the output and the figures.
- Following Example 1, solve the following IVP:
x
0 = 2x + 8y
y
0 = −x − 2y,
with the initial condition x(0) = 2, y(0) = −1. Plot the solutions and the phase portrait. - Following Example 2, find the numerical solution to Problem 2 above for 0 ≤ x ≤ 10 with increments h = 0.05
(you can try using even smaller increment). Plot the solutions and the phase portrait. - (optional with bonus) Compare Problems 2 and 3, why the phase portraits behave differently here?
5
MATH 430 Lab 2 Fall 2019 - Following Example 5 to solve the following linear system:
x
0 =
−2 0
0 −2
x +
2e
−2t
3t
, x(0) =
1
−1
,
• Save all your command inputs, command window output, and figures to a script file.
• The assignment need to be handed in a single file. Here is one way to do it. Open a blank Window Word
file. Copy and paste the required command inputs, their outputs from Matlab’s Command Window, together
with the figures. For the required figures, save them as JPEG files or any graphic files and then import/insert
them to the Word file. Organize your document under example/problem headings, Example 1, 2, 3, Problem
1, 2, 3, etc. Make sure to include your name, the lab assignment number, and class information. Save the
Word document as a PDF file. If you do not use Word, use any other program you prefer to accomplish the
same task.
• You may adjust the font size and spacing between paragraphs to save pages.
• Upload the scripts, results and figures in one single PDF file to Canvas.
6
MATH 430 Lab 2 Fall 2019
Expected figrues in your submission:
Figure 3: Plot of Example 1. Left are solutions. Right is the phase portrait.
Figure 4: Plot of Example 2. Left are solutions. Right is the phase portrait.
Figure 5: Plot of Problem 2. Left are solutions. Right is the phase portrait.
7
MATH 430 Lab 2 Fall 2019
Figure 6: Plot of Problem 3. Left are solutions. Right is the phase portrait.
8