n = 6; M = rand(n,n); M < .12
ans =
0 0 0 0 0 0
0 0 0 0 0 1
0 0 0 0 0 0
0 0 0 1 0 1
0 0 0 0 0 1
1 0 0 0 0 0
P = M < .12;
P * M
ans =
0 0 0 0 0 0
0.0975 0.9706 0.9157 0.9340 0.1712 0.8235
0 0 0 0 0 0
1.0109 1.9355 1.0576 0.9697 0.5634 0.8696
0.0975 0.9706 0.9157 0.9340 0.1712 0.8235
0.8147 0.2785 0.9572 0.7922 0.6787 0.7060
P .* M
ans =
0 0 0 0 0 0
0 0 0 0 0 0.0318
0 0 0 0 0 0
0 0 0 0.0357 0 0.0462
0 0 0 0 0 0.0971
0.0975 0 0 0 0 0
(P .* M) \ [1:6]'
[Warning: Matrix is singular to working precision.]
ans =
NaN
NaN
NaN
NaN
-Inf
Inf
sparse(P .* M)
ans =
(6,1) 0.0975
(4,4) 0.0357
(2,6) 0.0318
(4,6) 0.0462
(5,6) 0.0971
B = sparse(P .* M)
B =
(6,1) 0.0975
(4,4) 0.0357
(2,6) 0.0318
(4,6) 0.0462
(5,6) 0.0971
B
B =
(6,1) 0.0975
(4,4) 0.0357
(2,6) 0.0318
(4,6) 0.0462
(5,6) 0.0971
spy(B)
n = 100; M = rand(n,n); spy(M .* (M < .03))
close all; n = 1000; M = rand(n,n); spy(M .* (M < .03))
n = 6; M = rand(n,n); spy(M .* (M < .03))
n = 6; M = rand(n,n); B = M .* (M < .12);
B
B =
0 0.0867 0 0 0 0.0435
0 0 0 0.0788 0 0
0 0 0 0 0.0846 0
0 0 0 0 0 0.1055
0 0 0 0 0.0634 0.0933
0 0 0 0 0 0
B = sparse(B)
B =
(1,2) 0.0867
(2,4) 0.0788
(3,5) 0.0846
(5,5) 0.0634
(1,6) 0.0435
(4,6) 0.1055
(5,6) 0.0933
B + eye(n)
ans =
1.0000 0.0867 0 0 0 0.0435
0 1.0000 0 0.0788 0 0
0 0 1.0000 0 0.0846 0
0 0 0 1.0000 0 0.1055
0 0 0 0 1.0634 0.0933
0 0 0 0 0 1.0000
sparse(B + eye(n))
ans =
(1,1) 1.0000
(1,2) 0.0867
(2,2) 1.0000
(3,3) 1.0000
(2,4) 0.0788
(4,4) 1.0000
(3,5) 0.0846
(5,5) 1.0634
(1,6) 0.0435
(4,6) 0.1055
(5,6) 0.0933
(6,6) 1.0000
B + sparse(eye(n))
ans =
(1,1) 1.0000
(1,2) 0.0867
(2,2) 1.0000
(3,3) 1.0000
(2,4) 0.0788
(4,4) 1.0000
(3,5) 0.0846
(5,5) 1.0634
(1,6) 0.0435
(4,6) 0.1055
(5,6) 0.0933
(6,6) 1.0000
sparse(B + eye(n)) \ [1:n]'
ans =
0.5886
1.7348
2.6467
3.3671
4.1757
6.0000
(B + eye(n)) \ [1:n]'
ans =
0.5886
1.7348
2.6467
3.3671
4.1757
6.0000
sparse(B + eye(n))
ans =
(1,1) 1.0000
(1,2) 0.0867
(2,2) 1.0000
(3,3) 1.0000
(2,4) 0.0788
(4,4) 1.0000
(3,5) 0.0846
(5,5) 1.0634
(1,6) 0.0435
(4,6) 0.1055
(5,6) 0.0933
(6,6) 1.0000
fulla(sparse(B + eye(n)))
{Undefined function or variable 'fulla'.
}
full(sparse(B + eye(n)))
ans =
1.0000 0.0867 0 0 0 0.0435
0 1.0000 0 0.0788 0 0
0 0 1.0000 0 0.0846 0
0 0 0 1.0000 0 0.1055
0 0 0 0 1.0634 0.0933
0 0 0 0 0 1.0000
pwd
ans =
/Users/pauldj/works/courses/LSC/LSC-17/FPA
cd ../Matlab/
edit ode1
a=0; b=3; N=10; h=(b-a)/N
h =
0.3000
a=0; b=3; N=10; h=(b-a)/N;
t=linspace(a, b, N+1)
t =
Columns 1 through 7
0 0.3000 0.6000 0.9000 1.2000 1.5000 1.8000
Columns 8 through 11
2.1000 2.4000 2.7000 3.0000
u(1) = 1;
for i=1:N, ...
u(i+1) = u(i) + h * ode1( t(i), u(i) );
end
u
u =
Columns 1 through 7
1.0000 0.7000 -0.6184 -0.5491 0.2118 0.2745 -0.1218
Columns 8 through 11
-0.1874 0.0304 0.0943 -0.0150
t
t =
Columns 1 through 7
0 0.3000 0.6000 0.9000 1.2000 1.5000 1.8000
Columns 8 through 11
2.1000 2.4000 2.7000 3.0000
plot( t, u )
a=0; b=3; N=100; h=(b-a)/N; close all; t=linspace(a,b,N+1); u(1)=1;
for i=1:N, u(i+1) = u(i) + h * ode1( t(i), u(i) ); end
plot( t, u )
help ode45
ode45 Solve non-stiff differential equations, medium order method.
[TOUT,YOUT] = ode45(ODEFUN,TSPAN,Y0) with TSPAN = [T0 TFINAL] integrates
the system of differential equations y' = f(t,y) from time T0 to TFINAL
with initial conditions Y0. ODEFUN is a function handle. For a scalar T
and a vector Y, ODEFUN(T,Y) must return a column vector corresponding
to f(t,y). Each row in the solution array YOUT corresponds to a time
returned in the column vector TOUT. To obtain solutions at specific
times T0,T1,...,TFINAL (all increasing or all decreasing), use TSPAN =
[T0 T1 ... TFINAL].
[TOUT,YOUT] = ode45(ODEFUN,TSPAN,Y0,OPTIONS) solves as above with default
integration properties replaced by values in OPTIONS, an argument created
with the ODESET function. See ODESET for details. Commonly used options
are scalar relative error tolerance 'RelTol' (1e-3 by default) and vector
of absolute error tolerances 'AbsTol' (all components 1e-6 by default).
If certain components of the solution must be non-negative, use
ODESET to set the 'NonNegative' property to the indices of these
components.
ode45 can solve problems M(t,y)*y' = f(t,y) with mass matrix M that is
nonsingular. Use ODESET to set the 'Mass' property to a function handle
MASS if MASS(T,Y) returns the value of the mass matrix. If the mass matrix
is constant, the matrix can be used as the value of the 'Mass' option. If
the mass matrix does not depend on the state variable Y and the function
MASS is to be called with one input argument T, set 'MStateDependence' to
'none'. ODE15S and ODE23T can solve problems with singular mass matrices.
[TOUT,YOUT,TE,YE,IE] = ode45(ODEFUN,TSPAN,Y0,OPTIONS) with the 'Events'
property in OPTIONS set to a function handle EVENTS, solves as above
while also finding where functions of (T,Y), called event functions,
are zero. For each function you specify whether the integration is
to terminate at a zero and whether the direction of the zero crossing
matters. These are the three column vectors returned by EVENTS:
[VALUE,ISTERMINAL,DIRECTION] = EVENTS(T,Y). For the I-th event function:
VALUE(I) is the value of the function, ISTERMINAL(I)=1 if the integration
is to terminate at a zero of this event function and 0 otherwise.
DIRECTION(I)=0 if all zeros are to be computed (the default), +1 if only
zeros where the event function is increasing, and -1 if only zeros where
the event function is decreasing. Output TE is a column vector of times
at which events occur. Rows of YE are the corresponding solutions, and
indices in vector IE specify which event occurred.
SOL = ode45(ODEFUN,[T0 TFINAL],Y0...) returns a structure that can be
used with DEVAL to evaluate the solution or its first derivative at
any point between T0 and TFINAL. The steps chosen by ode45 are returned
in a row vector SOL.x. For each I, the column SOL.y(:,I) contains
the solution at SOL.x(I). If events were detected, SOL.xe is a row vector
of points at which events occurred. Columns of SOL.ye are the corresponding
solutions, and indices in vector SOL.ie specify which event occurred.
Example
[t,y]=ode45(@vdp1,[0 20],[2 0]);
plot(t,y(:,1));
solves the system y' = vdp1(t,y), using the default relative error
tolerance 1e-3 and the default absolute tolerance of 1e-6 for each
component, and plots the first component of the solution.
Class support for inputs TSPAN, Y0, and the result of ODEFUN(T,Y):
float: double, single
See also ode23, ode113, ode15s, ode23s, ode23t, ode23tb, ode15i,
odeset, odeplot, odephas2, odephas3, odeprint, deval,
odeexamples, rigidode, ballode, orbitode, function_handle.
Reference page in Help browser
doc ode45
ode45( @ode1, [0 3], 1 )
hold on
plot( t, u, 'r+' )
sol3 = e^-3 * cos( 15 )
{Undefined function or variable 'e'.
}
sol3 = E^-3 * cos( 15 )
{Undefined function or variable 'E'.
}
sol3 = exp^-3 * cos( 15 )
{Error using exp
Not enough input arguments.
}
sol3 = exp(-3) * cos( 15 )
sol3 =
-0.0378
format long; sol3
sol3 =
-0.037822634055742
a=0; b=3; N=10; h=(b-a)/N; close all; t=linspace(a,b,N+1); u(1)=1;
for i=1:N, u(i+1) = u(i) + h * ode1( t(i), u(i) ); end
abs( (u(end) - sol3)/sol3 )
ans =
0.001946729465011
a=0; b=3; N=100; h=(b-a)/N; close all; t=linspace(a,b,N+1); u(1)=1;
for i=1:N, u(i+1) = u(i) + h * ode1( t(i), u(i) ); end
abs( (u(end) - sol3)/sol3 )
ans =
0.001946729465011
a=0; b=3; N=10; h=(b-a)/N; close all; t=linspace(a,b,N+1); u(1)=1;
for i=1:N, u(i+1) = u(i) + h * ode1( t(i), u(i) ); end
abs( (u(N) - sol3)/sol3 )
ans =
3.492646298883530
a=0; b=3; N=100; h=(b-a)/N; close all; t=linspace(a,b,N+1); u(1)=1;
for i=1:N, u(i+1) = u(i) + h * ode1( t(i), u(i) ); end
abs( (u(N) - sol3)/sol3 )
ans =
0.125747733573885
a=0; b=3; N=10; h=(b-a)/N; close all; t=linspace(a,b,N+1); u(1)=1;
for i=1:N, u(i+1) = u(i) + h * ode1( t(i), u(i) ); end
abs( (u(N+1) - sol3)/sol3 )
ans =
0.602534223273937
a=0; b=3; N=100; h=(b-a)/N; close all; t=linspace(a,b,N+1); u(1)=1;
for i=1:N, u(i+1) = u(i) + h * ode1( t(i), u(i) ); end
abs( (u(N+1) - sol3)/sol3 )
ans =
0.001946729465011
a=0; b=3; N=1000; h=(b-a)/N; close all; t=linspace(a,b,N+1); u(1)=1;
for i=1:N, u(i+1) = u(i) + h * ode1( t(i), u(i) ); end
abs( (u(N+1) - sol3)/sol3 )
ans =
7.272437623209290e-04
a=0; b=3; N=10000; h=(b-a)/N; close all; t=linspace(a,b,N+1); u(1)=1;
for i=1:N, u(i+1) = u(i) + h * ode1( t(i), u(i) ); end
abs( (u(N+1) - sol3)/sol3 )
ans =
7.799831062937499e-05
a=0; b=3; N=100000; h=(b-a)/N; close all; t=linspace(a,b,N+1); u(1)=1;
for i=1:N, u(i+1) = u(i) + h * ode1( t(i), u(i) ); end
abs( (u(N+1) - sol3)/sol3 )
ans =
7.852521557373846e-06
a=0; b=3; N=1000000; h=(b-a)/N; close all; t=linspace(a,b,N+1); u(1)=1;
for i=1:N, u(i+1) = u(i) + h * ode1( t(i), u(i) ); end
abs( (u(N+1) - sol3)/sol3 )
ans =
7.857789108854612e-07
edit MYeuler
MYeuler( @ode1, 0, 3, 1, 100 )
{Error: File: MYeuler.m Line: 3 Column: 5
Expression or statement is incomplete or incorrect.
}
[t,u] = MYeuler( @ode1, 0, 3, 1, 100 );
{Error: File: MYeuler.m Line: 3 Column: 5
Expression or statement is incomplete or incorrect.
}
pwd
ans =
/Users/pauldj/works/courses/LSC/LSC-17/Matlab
sqrts( 1.0 )
7, 1, 1
8, 1, 1
9, 1, 1
10, 1, 1
11, 1, 1
12, 1, 1
13, 1, 1
14, 1, 1
15, 1, 1
sqrts( 1.0 )
7, 8, 9, 10, 11, 12, 13, 14, 15,
sqrts( 1.0 )
7, 8, 9, 10, 11, 12, 13, 14, 15, sqrts( 1.0 )
7, 1.000000000000000, 1.000000000000000
8, 1.000000000000000, 1.000000000000000
9, 1.000000000000000, 1.000000000000000
10, 1.000000000000000, 1.000000000000000
11, 1.000000000000000, 1.000000000000000
12, 1.000000000000000, 1.000000000000000
13, 1.000000000000000, 1.000000000000000
14, 1.000000000000000, 1.000000000000000
15, 1.000000000000000, 1.000000000000000
sqrts( 1.0, 25 )
7, 1.000000000000000, 1.000000000000000
8, 1.000000000000000, 1.000000000000000
9, 1.000000000000000, 1.000000000000000
10, 1.000000000000000, 1.000000000000000
11, 1.000000000000000, 1.000000000000000
12, 1.000000000000000, 1.000000000000000
13, 1.000000000000000, 1.000000000000000
14, 1.000000000000000, 1.000000000000000
15, 1.000000000000000, 1.000000000000000
16, 1.000000000000000, 1.000000000000000
17, 1.000000000000000, 1.000000000000000
18, 1.000000000000000, 1.000000000000000
19, 1.000000000000000, 1.000000000000000
20, 1.000000000000000, 1.000000000000000
21, 1.000000000000000, 1.000000000000000
22, 1.000000000000000, 1.000000000000000
23, 1.000000000000000, 1.000000000000000
24, 1.000000000000000, 1.000000000000000
25, 1.000000000000000, 1.000000000000000
0000
ans =
0
sqrts( 2.0, 15 )
7, 1.999999999999971, 2.000000000000000
8, 2.000000000000024, 2.000000000000000
9, 2.000000000000024, 2.000000000000000
10, 2.000000000000024, Inf
11, 2.000000000000024, Inf
12, 1.999999999999134, Inf
13, 1.999999999997329, Inf
14, 1.999999999997329, Inf
15, 1.999999999997329, Inf
sqrts( 2.0, 20 )
7, 1.999999999999971, 2.000000000000000
8, 2.000000000000024, 2.000000000000000
9, 2.000000000000024, 2.000000000000000
10, 2.000000000000024, Inf
11, 2.000000000000024, Inf
12, 1.999999999999134, Inf
13, 1.999999999997329, Inf
14, 1.999999999997329, Inf
15, 1.999999999997329, Inf
16, 2.000000000011775, Inf
17, 2.000000000040858, Inf
18, 2.000000000040858, Inf
19, 2.000000000157359, Inf
20, 2.000000000157359, Inf
sqrts( .5, 20 )
7, 0.499999999999999, 0.500000000000000
8, 0.499999999999999, 0.500000000000000
9, 0.499999999999999, 0.500000000000000
10, 0.499999999999999, 0.500000000000000
11, 0.499999999999999, 0.000000000000000
12, 0.500000000000114, 0.000000000000000
13, 0.500000000000114, 0.000000000000000
14, 0.500000000000571, 0.000000000000000
15, 0.500000000001479, 0.000000000000000
16, 0.500000000001479, 0.000000000000000
17, 0.500000000001479, 0.000000000000000
18, 0.500000000008752, 0.000000000000000
19, 0.500000000008752, 0.000000000000000
20, 0.500000000037885, 0.000000000000000
sqrts( .1, 20 )
7, 0.099999999999999, 0.100000000000000
8, 0.099999999999999, 0.100000000000000
9, 0.100000000000002, 0.000000000000000
10, 0.100000000000002, 0.000000000000000
11, 0.100000000000002, 0.000000000000000
12, 0.100000000000002, 0.000000000000000
13, 0.100000000000002, 0.000000000000000
14, 0.100000000000095, 0.000000000000000
15, 0.100000000000095, 0.000000000000000
16, 0.100000000000460, 0.000000000000000
17, 0.100000000000460, 0.000000000000000
18, 0.100000000000460, 0.000000000000000
19, 0.100000000003375, 0.000000000000000
20, 0.099999999997554, 0.000000000000000
sqrts( .9, 20 )
7, 0.899999999999999, 0.900000000000000
8, 0.899999999999999, 0.900000000000000
9, 0.899999999999999, 0.900000000000000
10, 0.899999999999999, 0.900000000000000
11, 0.899999999999999, 0.900000000000000
12, 0.900000000000203, 0.900000000000000
13, 0.899999999999792, 0.000000000000000
14, 0.900000000000613, 0.000000000000000
15, 0.900000000000613, 0.000000000000000
16, 0.900000000000613, 0.000000000000000
17, 0.899999999994060, 0.000000000000000
18, 0.899999999994060, 0.000000000000000
19, 0.900000000020266, 0.000000000000000
20, 0.900000000072656, 0.000000000000000
sqrts( 4.1, 20 )
7, 4.099999999999977, 4.100000000000000
8, 4.100000000000105, 4.100000000000000
9, 4.100000000000337, Inf
10, 4.099999999999857, Inf
11, 4.100000000000778, Inf
12, 4.100000000000778, Inf
13, 4.100000000000778, Inf
14, 4.100000000008135, Inf
15, 4.099999999993023, Inf
16, 4.099999999993023, Inf
17, 4.100000000053106, Inf
18, 4.099999999933824, Inf
19, 4.100000000172273, Inf
20, 4.100000000648643, Inf
pwd
ans =
/Users/pauldj/works/courses/LSC/LSC-17/Matlab
ls
LQ.m diary-2.txt my_pi_up.m sqrtsq.m
MYeuler.m diary-3.txt myfilter.m sqsqrt.m
MYeuler.m~ hw13.m ode1.m traverse.m
Mv.m least.m scripts.tgz
diary-1.txt my_pi_down.m sqrts.m
edit MYlu
MYlu( 0.1 );
MYlu( 0.1 )
ans =
1 0
10 1
[L,U] = MYlu( 0.1 )
L =
1 0
10 1
U =
0.1000000 -1.0000000
0 11.0000000
[L,U] = MYlu( 0.1 );
eps = 0.1; [L,U] = MYlu( eps ); norm(L*U - [eps -1; 1 1])
ans =
1.4901161e-09
eps = 0.1; [L,U] = MYlu( eps ); norm(L*U - [eps -1; 1 1])/norm([eps -1; 1 1])
ans =
9.3643127e-10
eps = 0.1; [L,U] = MYlu( eps ); norm(L*U - [eps -1; 1 1])
ans =
1.4901161e-09
eps = 10^-4; [L,U] = MYlu( eps ); norm(L*U - [eps -1; 1 1])
ans =
2.5262125e-12
eps = 10^-8; [L,U] = MYlu( eps ); norm(L*U - [eps -1; 1 1])
ans =
1
eps = 10^-6; [L,U] = MYlu( eps ); norm(L*U - [eps -1; 1 1])
ans =
2.5247572e-15
eps = 10^-7; [L,U] = MYlu( eps ); norm(L*U - [eps -1; 1 1])
ans =
1.1686097e-15
eps = 10^-8; [L,U] = MYlu( eps ); norm(L*U - [eps -1; 1 1])
ans =
1
diary