Bài giảng Matlab ứng dụng - Phương trình vi phân thường

1.0
1.5
1.8
2.0
3.0
3.5
4.5
α maxρ
a) Luật Parabol
b) Luật tổ hợp tuyến tính 
32
2
1 cxcxc ++
xc
x
c
2
1 +
TS
N
g
u
y
ễ
n
H
o
à
i
S
ơ
n
Matlab program clear all
clc
alpha=[1 1.5 1.8 2.0 3.0 3.5 4.5]';
rho= [0.098158 0.075798 0.066604 0.049851
0.046624 0.04189 0.0346]';
% luật parabol qua 7 điểm: c1x^2+c2x+c3
A=[alpha.^2 alpha ones(size(alpha))];
disp(A'*A)
disp(A'*rho)
c =(A'*A)\(A'*rho)
% vẽ đồ thị
xfit=linspace(min(alpha),max(alpha));
yfit1=c(1)*xfit.^2+c(2)*xfit+c(3);
% luật c1/x+c2x
A=[1./alpha alpha];
c=(A'*A)\(A'*rho);
xfit=linspace(min(alpha),max(alpha));
yfit2=c(1)./xfit+c(2)*xfit;
plot(alpha,rho,'o',xfit,yfit1,'r',xfit,yfit2,'c')
xlabel('alpha')
ylabel('rho')
title(‘rho=f(alpha)')
legend(‘ dữ liệu đo đạc','luật parabol',luật tổ hợp')
grid on
TS
N
g
u
y
ễ
n
H
o
à
i
S
ơ
n
• II. Dùng tích phân số :
• 1. Luật hình thang (Trapzoidal Rule) :
))()((
2
)( 1
1
ii
x
x
i xfxf
h
dxxf
i
i
+≈ −∫
−
x0 = a x1 x2 …. Xn-1 xn=b
y
f(x)
x
Eh
E
xfxf
xfxfxfhI
nn
trap +⎟⎟⎠
⎞
⎜⎜⎝
⎛
++
+++=
− )()(2
.....)(2)(2)(
2 1
210
bxaxhiax
N
abh ni ==+=−= ,,*, 0
( ) EfffffhI nntrap ++++++= −1210 2.....222
Ví dụ
∑
=
+− +=−−≈
N
i
ii
ii
xxxxf
N
abE
1
11''
3
3
2
),()(
12
1
∫∫ ⎟⎟⎠
⎞
⎜⎜⎝
⎛ ⎟⎠
⎞⎜⎝
⎛+==
2
0
222
0 2
1)( dxxdxxfS πTính tích phân:
TS
N
g
u
y
ễ
n
H
o
à
i
S
ơ
n
Matlab program
clear all
clc
N=16;
a=0;
b=2;
h=(b-a)/N;
S=0;
for i=0:N
x=a+i*h;
if i==0 | i==N
c=1;
else
c=2;
end
S=S+c*pi*(1+(x/2).^2).^2;
end
S=h*S/2
Kết qủa:
-1.0341
-0.2609
-0.0654
-0.0163
-0.0040
-0.0010
12.7627
11.9895
11.7940
11.7449
11.7326
11.7296
1.
0.5
0.25
0.125
0.0625
0.03125
2
4
8
16
32
64
EhShhN
• 2. Luật Simpson 1/3 (Simpson Rule) :
[ ]∫ +++== b
a
EbfxfafhdxxfS )()(4)(
3
)(
2
,
2
,, 20
baxabhbxax +=−===
[ ]∫ +++== b
a
EfffhdxxfS 210 43
)(
EbfihafihafafhS
N
i
N
i
simp +⎥⎦
⎤⎢⎣
⎡ +++++= ∑ ∑−
=
−
=
1
1
2
2
)()(2)(4)(
3
E
xfxf
xfxfxfxfhS
nn
simp +⎟⎟⎠
⎞
⎜⎜⎝
⎛
++
++++=
− )()(4.....
)(4)(2)(4)(
3 1
3210
∑
=
+− +==−≈
N
i
ii
ii
xxxNxfffhNE
1
11''''''''''''
5
2
,/)(,
902
TS
N
g
u
y
ễ
n
H
o
à
i
S
ơ
n
Matlab program
clear all
clc
N=16;
a=0;
b=2;
h=(b-a)/N;
S=0;
for i=0:N
x=a+i*h;
if i==0 | i==N
c=1;
elseif i==fix(i/2)*2+1
c=4;
else
c=2;
end
S=S+c*pi*(1+(x/2).^2).^2;
end
S=h*S/3
-0.0523
-0.0032
-0.0002
-0.0000
-0.0000
-0.0000
11.7809
11.7318
11.7288
11.7286
11.7286
11.7286
1.
0.5
0.25
0.125
0.0625
0.03125
2
4
8
16
32
64
EhShhN
Kết qủa:
1.32
5
0 10 20 30 40 50
Số phân đoạn
1.32
1.31
5
1.31
1.30
5
1.3
1.29
5
1.29
Luật Simpson
Luật hình thang
Chính xác
Giá trị tích phân
TS
N
g
u
y
ễ
n
H
o
à
i
S
ơ
n
Sai số phương pháp:
TS
N
g
u
y
ễ
n
H
o
à
i
S
ơ
n
3. Tích phân Gauss (Gauss quadrature):
)()()()( 2211
1
1 nn
xfwxfwxfwdxxfI +++≈= ∫− L
clear all
clc
format long
x1=-0.861136;
x2=-0.339981;
x3=0.339981;
x4=0.861136;
% ------trọng số------------------------------------
w1=0.347855;
w2=0.652145;
w3=0.652145;
w4=0.347855;
f1=w1*gauss1(x1);
f2=w2*gauss1(x2);
f3=w3*gauss1(x3);
f4=w4*gauss1(x4);
m=f1+f2+f3+f4
%--------------------------------------------------------------------
function ff=gauss1(x)
ff=400*x^5-900*x^4+675*x^3-200*x^2+25*x+0.2;
kết quả:
I=-4.929329328775451e+002
Ví dụ: ∫
−
+−+−+=
1
1
5432 )400900675200252.0( dxxxxxxI
Tính với 4 điểm Gauss:
Matlab program
TS
N
g
u
y
ễ
n
H
o
à
i
S
ơ
n
MATLAB - FEM
Bài tập 3.4
clear all; clc; close all
echo off
%-------------------------------------------------------------
Edof=[1 1 2 3 4 5 6;
2 4 5 6 7 8 9];
%-------------------------------------------------------------
K=zeros(9); f=zeros(9,1); f(8)=-88.9/2;
%-------------------------------------------------------------
h=17.9; tw=0.315; bf=6.015;tf=0.525;
A=2*tf*bf+tw*(h-2*tf);
I=2.5e-2; E=2.1e8; L=6.1;
ep=[E A I];
Ex=[0 L;
L 3*L/2];
Ey=zeros(2,2);
Eq=zeros(2,2);
%-------------------------------------------------------------
for i=1:2
[Ke,fe]=beam2e(Ex(i,:),Ey(i,:),ep);
[K,f]=assem(Edof(i,:),K,Ke,f,fe);
end
Ed=extract(Edof,a);
[es1,edi1,eci1]=beam2s(Ex(1,:),Ey(1,:),ep,Ed(1,:),Eq(1,:),20);
[es2,edi2,eci2]=beam2s(Ex(2,:),Ey(2,:),ep,Ed(2,:),Eq(2,:),10);
%-------------------------------------------------------------
%-------------------------------------------------------------
bc=[1 0;2 0;4 0;5 0;7 0;9 0];
a=solveq(K,f,bc);
%-------------------------------------------------------------
function [Ke,fe]=beam2e(ex,ey,ep,eq);
%---------------------------------------------------------------------
% INPUT: 
% ex = [x1 x2]
% ey = [y1 y2] element node coordinates
%
% ep = [E A I] element properties
% E: Young's modulus
% A: Cross section area
% I: Moment of inertia
%
% eq = [qx qy] distributed loads, local directions
% 
% OUTPUT: Ke : element stiffness matrix (6 x 6)
% fe : element load vector (6 x 1)
%-------------------------------------------------------------
b=[ ex(2)-ex(1); ey(2)-ey(1) ];
L=sqrt(b'*b); n=b/L;
TS
N
g
u
y
ễ
n
H
o
à
i
S
ơ
n
MATLAB - FEM
E=ep(1); A=ep(2); I=ep(3);
qx=0; qy=0; if nargin>3; qx=eq(1); qy=eq(2); end
Kle=[E*A/L 0 0 -E*A/L 0 0 ;
0 12*E*I/L^3 6*E*I/L^2 0 -12*E*I/L^3 
6*E*I/L^2;
0 6*E*I/L^2 4*E*I/L 0 -6*E*I/L^2 2*E*I/L;
-E*A/L 0 0 E*A/L 0 0 ;
0 -12*E*I/L^3 -6*E*I/L^2 0 12*E*I/L^3 -
6*E*I/L^2;
0 6*E*I/L^2 2*E*I/L 0 -6*E*I/L^2 
4*E*I/L];
fle=L*[qx/2 qy/2 qy*L/12 qx/2 qy/2 -qy*L/12]';
G=[n(1) n(2) 0 0 0 0;
-n(2) n(1) 0 0 0 0;
0 0 1 0 0 0;
0 0 0 n(1) n(2) 0;
0 0 0 -n(2) n(1) 0;
0 0 0 0 0 1];
Ke=G'*Kle*G; fe=G'*fle; 
%--------------------------end------------------------------
function [K,f]=assem(edof,K,Ke,f,fe)
%-------------------------------------------------------------
% INPUT: edof: dof topology matrix
% K : the global stiffness matrix
% Ke: element stiffness matrix
% f : the global force vector
% fe: element force vector
%
% OUTPUT: K : the new global stiffness matrix
% f : the new global force vector
%-------------------------------------------------------------
[nie,n]=size(edof);
t=edof(:,2:n);
for i = 1:nie
K(t(i,:),t(i,:)) = K(t(i,:),t(i,:))+Ke;
if nargin==5
f(t(i,:))=f(t(i,:))+fe;
end
end
%--------------------------end--------------------------------
TS
N
g
u
y
ễ
n
H
o
à
i
S
ơ
n
MATLAB - FEM
function [d,Q]=solveq(K,f,bc)
% a=solveq(K,f)
% [a,Q]=solveq(K,f,bc)
%-------------------------------------------------------------
% PURPOSE
% Solve static FE-equations considering boundary conditions.
%
% INPUT: K : global stiffness matrix, dim(K)= nd x nd
% f : global load vector, dim(f)= nd x 1
%
% bc : boundary condition matrix
% dim(bc)= nbc x 2, nbc : number of b.c.'s
%
% OUTPUT: a : solution including boundary values
% Q : reaction force vector
% dim(a)=dim(Q)= nd x 1, nd : number of dof's
%-------------------------------------------------------------
if nargin==2 ; 
d=K\f ; 
elseif nargin==3;
[nd,nd]=size(K);
fdof=[1:nd]';
%
d=zeros(size(fdof));
Q=zeros(size(fdof));
%
pdof=bc(:,1);
dp=bc(:,2);
fdof(pdof)=[];
s=K(fdof,fdof)\(f(fdof)-K(fdof,pdof)*dp);
%A=K(fdof,fdof);
%B=(f(fdof)-K(fdof,pdof)*dp);
%s=pcg(A,B);
%
d(pdof)=dp;
d(fdof)=s;
end 
Q=K*d-f;
%--------------------------end--------------------------------
function [ed]=extract(edof,a)
%-------------------------------------------------------------
% PURPOSE
% Extract element displacements from the global 
% displacement
% vector according to the topology matrix edof.
% INPUT: a: the global displacement vector
% edof: topology matrix
% OUTPUT: ed: element displacement matrix
%-------------------------------------------------------------
[nie,n]=size(edof);
t=edof(:,2:n);
%
for i = 1:nie
ed(i,1:(n-1))=a(t(i,:))';
end
%
%--------------------------end--------------------------------
TS
N
g
u
y
ễ
n
H
o
à
i
S
ơ
n
MATLAB - FEM
function [es,edi,eci]=beam2s(ex,ey,ep,ed,eq,n)
% PURPOSE
% Compute section forces in two dimensional beam element 
%(beam2e). 
% INPUT: ex = [x1 x2]
% ey = [y1 y2] element node coordinates
% ep = [E A I] element properties,
% E: Young's modulus
% A: cross section area
% I: moment of inertia
% ed = [u1 ... u6] element displacements
% eq = [qx qy] distributed loads, local directions 
% n : number of evaluation points ( default=2 )
% 
% OUTPUT: 
% es = [ N1 V1 M1 ; section forces, local directions, in N2 V2 M2 
%; n points along the beam, dim(es)= n x 3 .........] 
% edi = [ u1 v1 ; element displacements, local directions, u2 v2 
%; in n points along the beam, dim(es)= n x 2 ....] 
% eci = [ x1 ; local x-coordinates of the evaluation 
% x2 ; points, (x1=0 and xn=L) ...]
%-------------------------------------------------------------
EA=ep(1)*ep(2); EI=ep(1)*ep(3);
b=[ ex(2)-ex(1); ey(2)-ey(1) ];
L=sqrt(b'*b);
C=[0 0 0 1 0 0;
0 0 0 0 0 1;
0 0 0 0 1 0;
L 0 0 1 0 0;
0 L^3 L^2 0 L 1;
0 3*L^2 2*L 0 1 0];
n=b/L;
G=[n(1) n(2) 0 0 0 0;
-n(2) n(1) 0 0 0 0;
0 0 1 0 0 0;
0 0 0 n(1) n(2) 0;
0 0 0 -n(2) n(1) 0;
0 0 0 0 0 1];
M=inv(C)*(G*ed'-[0 0 0 -qx*L^2/(2*EA) 
qy*L^4/(24*EI) qy*L^3/(6*EI)]' );
A=[M(1) M(4)]'; B=[M(2) M(3) M(5) M(6)]';
x=[0:L/(ne-1):L]'; zero=zeros(size(x)); 
one=ones(size(x));
u=[x one]*A-(x.^2)*qx/(2*EA);
du=[one zero]*A-x*qx/EA;
v=[x.^3 x.^2 x one]*B+(x.^4)*qy/(24*EI); 
% dv=[3*x.^2 2*x one zero]*B+(x.^3)*qy/(6*EI);
d2v=[6*x 2*one zero zero]*B+(x.^2)*qy/(2*EI);
d3v=[6*one zero zero zero]*B+x*qy/EI;
N=EA*du; M=EI*d2v; V=-EI*d3v; 
es=[N V M];
edi=[u v];
eci=x;
%--------------------------end--------------------------------
if length(ed(:,1)) > 1 
disp('Only one row is allowed in the ed matrix !!!')
return
end
qx=0; qy=0; if nargin>4; qx=eq(1); qy=eq(2); end 
ne=2; if nargin>5; ne=n; end;