Download MENG274 Numerical Codes and more Exercises Numerical Methods in Engineering in PDF only on Docsity!
GUASS ELIMINATION
A=[2 -3 -1;3 2 -5;2 4 -1];
b=[3 -9 -5]'; Ab=[A b]; n=length(A); i=2;k=1;Ab(i,:)=Ab(i,:)-(Ab(i,k)/Ab(k,k))Ab(k,:); i=3;k=1;Ab(i,:)=Ab(i,:)-(Ab(i,k)/Ab(k,k))Ab(k,:); i=3;k=2;Ab(i,:)=Ab(i,:)-(Ab(i,k)/Ab(k,k))*Ab(k,:); U=Ab(:,1:3); c=Ab(:,4); x=U\c A\b %swap rows using A=[1 3 4; 2 4 5]; A([1 2],:)=A([2 1],:)
LU DECOMPOSITION
A=[-3 6 -4;9 -8 24;-12 24 -26];
b=[-3 65 -42]'; n=length(A); L=eye(n); U=A; i=2;k=1; L(i,k)=A(i,k)/A(k,k);U(i,:)=U(i,:)-L(i,k)U(k,:); i=3;k=1; L(i,k)=A(i,k)/A(k,k);U(i,:)=U(i,:)-L(i,k)U(k,:); i=3;k=2; L(i,k)=A(i,k)/A(k,k);U(i,:)=U(i,:)-L(i,k)*U(k,:); y=L\b; x=U\y A\b
CHOLESKI
A=[1 1 1;1 2 2;1 2 3];
b=[1 3/2 3]'; L=[whatever you solved by hand]; y=L\b; x=L'\y %IN COMMAND WINDOW syms L11 L21 L31 L22 L32 L33 real L=[L11 0 0;L21 L22 0;L31 L32 L33] L = [ L11, 0, 0] [ L21, L22, 0] [ L31, L32, L33] LL' ans = [ L11^2, L11L21, L11L31] [ L11L21, L21^2 + L22^2, L21L31 + L22L32]
[ L11L31, L21L31 + L22*L32, L31^2 + L32^2 + L33^2]
CURVE FITTING(QUADRATIC PLUS LINEAR)
x=[1 2.5 3.5 4 1.1 1.8 2.2 3.7]';
y=[6.008 15.722 27.120 33.772 5.257 9.549 11.098 28.828]';
n=length(x);
X_L=[ones(n,1) x];
c_L=inv(X_L'X_L)X_L'*y;
m_L=length(c_L);
a_L=c_L(1);b_L=c_L(2);
f_L=@(x) a_L+b_L*x;
e1=y-f_L(x);
s_1=sum(e1.^2);
sigma_L=sqrt(s_1/(n-m_L))
X_Q=[ones(n,1) x x.^2];
c_Q=inv(X_Q'X_Q)X_Q'*y;
m_Q=length(c_Q);
a_Q=c_Q(1);b_Q=c_Q(2);c=c_Q(3);
f_Q=@(x) a_Q+b_Qx+cx.^2;
e2=y-f_Q(x);
s=sum(e2.^2);
sigma_Q=sqrt(s/(n-m_Q))
hold all
fplot(f_L)
fplot(f_Q)
plot(x,y,'ro')
BISECTION
syms x
ezplot(function)
grid on
a=0;b=1;
f=@(x) function
[a b f(a) f((a+b)/2) f(b)]
NEWTON RAPHSON (NON-LINEAR SYSTEMS)
f=sin(x)+3*cos(x)-2;
df=cos(x) - 3*sin(x);
x=x-f/df;
i=i+
syms x y
sol=dsolve('Dy+4*y==x^2','y(0)==1','x');
subs(sol,x,0.1);
vpa(ans)
IVP(EULER)
x=0;y=1;h=0.1;
f=@(x,y) x^2-4*y;
xf=10;
xSol=x;
ySol=y;
i=1;
while x<=xf
xSol(i)=x;
ySol(i)=y;
y=y+h*f(x,y);
x=x+h;
i=i+1;
end
[xSol' ySol']
IVP(RK2)
t=0;y=1;h=0.1;
f=@(t,y) t^2-4*y;
tf=0.1;
i=1;
while t<=tf
tSol(i)=t;
ySol(i)=y;
k1=h*f(t,y);
k2=h*f(t+(h/2),y+(k1/2));
y=y+k2;
t=t+h;
i=i+1;
end
[tSol' ySol']
IVP(RK4 ORDER 1)
x=0;y=1;h=0.1;
f=@(x,y) x^2-4*y;
xf=0.1;
i=1;
while x<=xf
ysol(i)=y;
k1=h*f(x,y);
k2=h*f(x+(h/2),y+(k1/2));
k3=h*f(x+(h/2),y+(k2/2));
k4=h*f(x+h,y+k3);
y=y+1/6(k1+2k2+2*k3+k4);
x=x+h;
i=i+1;
end
[ysol']
IVP(RK4 ORDER 2 OR MORE)
x=0;y(1)=0;y(2)=1;h=0.25; f=@(x,y) [y(2) -0.1y(2)-x]; xf=2; i=1; while x<=xf ysol(i,:)=y; xsol(i,:)=x; k1=hf(x,y); k2=hf(x+(h/2),y+(k1/2)); k3=hf(x+(h/2),y+(k2/2)); k4=hf(x+h,y+k3); y=y+(1/6)(k1+2k2+2k3+k4); x=x+h; i=i+1; end [xsol ysol] plot(xsol,ysol)
BVP(DIRICHLET)
x=linspace(0,pi/2,n); h=x(2)-x(1); A=zeros(n); b=zeros(n,1); for i=2:n- A(i,i)=4h^2-2; A(i,i-1)=1; A(i,i+1)=1; b(i)=4h^2*x(i); end A(1,1)=1; b(1)=0; A(n,n)=1; b(n)=0; y=A\b plot(x,y)
BVP(NEUMAN)
x=linspace(0,pi/2,n); h=x(2)-x(1); A=zeros(n);
