Cheatography

# Numerical Methods 262 Cheat Sheet (DRAFT) by Kamva

Doing numerical methods using MATLAB to compute the certain codes needed.

This is a draft cheat sheet. It is a work in progress and is not finished yet.

### Curve fitting: Higher order Polyno­mials

The formula: ax + by + cz = 1
``A= [-1 -1 -3;1 0 0;2 1 2;2 3 6];``

``b = [1;1;1;1];``

``x = A\b;``

NOTE: Ax = b...backsl­ash(\) to solve for a,b,c

### Curve fitting: linear­ization

Expone­ntial linearized function
The formula is y = bemx:
``x = [1 2 3 4]; y = [6 3 2 1];   ``
Given
``a = polyfi­t(x­,lo­g(y),1)``

``>>a = -0.5781 2.3411; ``

``m = a(1) ``

``>>m = -0.5781;``

``b = exp(a(2));``

``>>b = 10.3923;``

The final product: y = 10.39e-0.5781x

### Eigenv­alues & Eigenv­ectors

``A=[2 -1 0;-1 2 -1;0 -1 2]; a=eig(A);``

``[V,D] = eig(A)``
V - eigenv­ectors
``V(:,1)``
:get 1st column
``v1 = V(:,1);``

``v1 = v1/v1(1);``
: normalise eigenv­ector

### ODEs using ode45

``a = 0.1; b = 0.002; c = 0.0025; d = 0.2;``

``f = @(t,x) [-ax(1)+bx(1)*x(2);``

``dx(2)-cx(1)*x­(2)];``

``X0 = [20; 100];``

``Tspan = linspa­ce(­0,1­50,­1000);``

``[T,X] = ode45(­f,T­spa­n,X0);``

``plot(T­,X,­'Li­neW­idt­h',2);``

### ODEs using ode45

``a = 0.1; b = 0.002; c = 0.0025; d = 0.2;``

``f = @(t,x) [-ax(1)+bx(1)*x(2);``

``dx(2)-cx(1)*x­(2)];``

``X0 = [20; 100];``

``Tspan = linspa­ce(­0,1­50,­1000);``

``[T,X] = ode45(­f,T­spa­n,X0);``

``plot(T­,X,­'Li­neW­idt­h',2);``

### Gauss Integr­ation

 Gauss 2-pt rule Gauss 3-pt rule
``f=@(x) 1/2*sq­rt(­16+­(x+­3).^­4)./(­x+3)^2;``

``G2=f(-­1/s­qrt­(3)­)+f­(1/­sqr­t(3));``

``>>G2= 1.1289;``

and
``G3=5/9f(-sqr­t(3­/5)­)+8/9f(0)+5­/9*­f(s­qrt­(3/5));``

``>>G3= 1.1319;``

### ODEs: Bounda­ry-­value problems (BVPs)

 Finite Difference Method
``y0=1;y1=2; % N = 1/h;h = 0.2``

``i = [1:N-1]';``

``xi = i*h;    % Interior nodes``

``A=(h^2­-2)­*ey­e(N-1); % Diagonal part of A ``

``A=A+di­ag(­1-0.5hxi(2:N­-1)­,-1);% Add sub- and superd­iagonal``

``A=A+di­ag(­1+0.5hxi(1:N­-2)­,+1);% Add sub- and superd­iagonal``

``b = 2h^2xi;% RHS refined formula``

``b(1)=b­(1)­-(1-0.5hxi(1))­*y0;% Add boundary contri­butions``

``b(N-1)­=b(­N-1­)-(­1+0.5hxi(N-1­))*y1;% Add boundary contri­butions``

``y = sparse­(A)\b;   % Solve``

### Finite Differ­ence: Laplac­e&­Poisson PDEs

You then use backslash after you have the equations

### Multiple Integrals

``f = @(x,y) x.^2-y.^2-1;``

``a = 0; b = 3; g = 0; p = 1;``

``integr­al2­(f,­a,b­,g,p);``

``>>ans = 5.0000­000­000­01309;``

or
``f = @(x,y) y;``

``a = 0; b = pi/4;``

``g = @(x) sin(x); p = @(x) cos(x);``

``m = integr­al2­(f,­a,b­,g,p);``

``>>m = 0.2500­000­000­00013;``

### Curve fittin­g& Interp­olation

 Least Squares: Backslash & Polyfit
y = a0 + a1x + a2x2:
``x = [1; 2; 3; 4]; y = [6; 3; 2; 1];``

``A = [ones(4,1) x x.ˆ2]; a = A\y``

``a = 9.5000 -4.1000 0.5000``

``a2 = polyfi­t(x­,y,2)``

``xx = linspa­ce(­0,5);``

``y2 = polyva­l(a­2,xx);``

``plot(x­,y,­'r.'­,x­x,y­2,'b-')``

### Error,­Res­idu­al,­Nor­ms,­Con­dition number

``x = [1; -1; 3; -5];``

``norm(x­,inf)``
||x||i­nf(­sub­script)
||A||2 = (||Ax||2)/(||x||2)
~ cond(A) = (||A||­)(||A-1||)
~
``cond(A)``
>> 1:ill-­con­dit­ion­ed&close singular matrix
~
``cond(A)``
≈ 1:well­-co­ndi­tioned
~ decimals digits that trustw­orthy = 10-a = 10c x 10-d
a = d-c
~r = b - Ax
det (A) = 0; No A-1; Singular matrix
test singul­arity:
``cond,c­ond­est­,rank``

### Built-in Functions

 Integral Trapz
``f=@(x)­sqr­t(1­+x.^­4)./x.^2;``

``L=inte­gra­l(f­,1,2);``

``>>L­=1.1­32­090­393­305918;``

or
``y=[0 0.3 0.5 1 1.5 2 2.5 3 4 5];``

``v=[0 0.4 0.5 0.56 0.6 0.63 0.66 0.68 0.71 0.72];``

``Q=10*t­rap­z(y,v);``

``>>Q­=30.8000;``

### Euler's Method

``f=@(x,y) -x*y;                     % Define rhs``

``N=100;                               % Number of steps``

``a=1; b=2;                           % Interval``

``h=(b-a)/N;                           % step size``

``x(1)=1; y(1)=4;                   % Initial values``

``for i=1:N``

``y(i+1)­=y(­i)+­h*f­(x(­i),­y(i));     % Euler``

``x(i+1)­=x(­i)+h;``

``end``

### Runge-­Kutta Methods

``f = @(x,y) -x*y;                   %Define rhs``

``N = 100;                             % Number of steps ``

....
``for i = 1:N``

``K1 = f(x(i)­,y(i));``

``K2 = f(x(i)+0.5h,y(i)+0.5h*K1);``

``K3 = f(x(i)+0.5h,y(i)+0.5h*K2);``

``K4 = f(x(i+­1),­y(i­)+h­*K3);``

``x(i+1) = x(i)+h;``

``y(i+1) = y(i)+(1/6)h(K1+2K­2+2­K3+K4);``

``end``

### Runge-­Kutta Methods

``f = @(x,y) -x*y;                   %Define rhs``

``N = 100;                             % Number of steps ``

....
``for i = 1:N``

``K1 = f(x(i)­,y(i));``

``K2 = f(x(i)+0.5h,y(i)+0.5h*K1);``

``K3 = f(x(i)+0.5h,y(i)+0.5h*K2);``

``K4 = f(x(i+­1),­y(i­)+h­*K3);``

``x(i+1) = x(i)+h;``

``y(i+1) = y(i)+(1/6)h(K1+2K­2+2­K3+K4);``

``end``

### Numerical Differ­ent­iation

``x = 0;``

``f = @(x) sqrt(1­-2*­sin­(x));``

``h = [0.002 0.001 0.0005 0.00025];``

``D2f = @(h) (f(x)-2f(x+h)­+f(x+2h))./h.^2;``

``D2f(h) = -1.0040 -1.0020 -1.0010 -1.0005; ``
*error decreases by factor 2 thus n = 1

### Composite Simpson's Rule

 Remember: N even
``f = @(x) sqrt(1­+x.^­4)./x.^2;``

``a = 1; b = 2;``

``N = 100; h = (b-a)/N;``

``x = a+[0:N]*h;``

``fx = f(x);``

``L = h/3(fx(1)+4sum(fx­(2:­2:N­))+­2*s­um(­fx(­3:2­:N-­1))­+fx­(N+1))``

or
``w = ones(s­ize­(x)); w(2:2:N) = 4; w(3:2:N-1) = 2;``

``L = h/3sum(w.f(x))``

``>>L = 1.1320­903­946­95845;``

### Composite Trapezium Rule

``f = @(x) sqrt(1­+x.^­4)./x.^2;``

``a = 1; b = 2;``

``N = 100; h = (b-a)/N;``

``x = a+[0:N]*h;``

``fx = f(x);``

``L = h3( 0.5fx(1) + sum(fx­(2:N)) + 0.5*fx­(N+1) )``

or
``w = ones(s­ize­(x)); w(1) = 0.5; w(N+1) = 0.5;``

``L = hsum(w.f(x));``

``>>L = 1.1321­016­727­88808;``

### Numerical Partial Differ­ent­iation

``x = 0; y = 0;``

``f = @(x,y) sqrt(1+2x-3y.^2);``

``h = [0.2 0.1 0.05 0.025];``

``Lapf = @(h) (f(x+h­,y)­+f(­x-h­,y)­+f(­x,y­+h)­+f(­x,y­-h)­-4*­f(x­,y)­)./­h.^2; ``

``Lapf(h) = -4.1505 -4.0356 -4.0088 -4.0022;``

### Modified Euler

``f = @(x,y) -x*y;    % Define rhs``

``N = 100;           % Number of steps``

``a = 1; b = 2;      % Interval``

``h= (b-a)/N;         % step size``

``x(1)=1­;y(1) = 4;  % Initial values``

``for i = 1:N7``

``yeu = y(i)+h­*f(­x(i­),y­(i)); % Euler``

``x(i+1) = x(i)+h;``

``y(i+1) = y(i)+0.5h(f(x(i­),y­(i)­)+f­(x(­i+1­),y­eu));``

``end``

### Spline

 Extrap­olate data set
``x = [22 35 48 61 74];``

``y = [320 490 540 500 480];``

``yy = spline­(x,­y,42);``
not-a-knot
``>>yy =  531.1971;``

``xx = linspa­ce(­20,­80,­1000);``

*clamped condit­ions:
``[0 y 0] ``

``ss = spline­(x,­y,xx);``

``plot(x­,y,­'r.'­,x­x,s­s,'­b-');``