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