-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathEx1292.m
144 lines (109 loc) · 4.72 KB
/
Ex1292.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
%----------------------------------------------------------------------------
% Example 12.9.2
% to solve the ordinary differential equation given as
% u'' - alpa*u*u'= 0, 0 < x < 1
% u(0) = 0 and u(1) = 1
% using quasilinearization.
%
% Variable descriptions
% k = element matrix
% f = element vector
% kk = system matrix
% ff = system vector
% index = a vector containing system dofs associated with each element
% bcdof = a vector containing dofs associated with boundary conditions
% bcval = a vector containing boundary condition values associated with
% the dofs in 'bcdof'
%----------------------------------------------------------------------------
%------------------------------------
% input data for control parameters
%------------------------------------
clear
nel=100; % number of elements
nnel=2; % number of nodes per element
ndof=1; % number of dofs per node
nnode=nel+1; % total number of nodes in system
sdof=nnode*ndof; % total system dofs
alpa=100; % coefficient of the nonlinear term
toler=0.0001; % error tolerance to terminate iterations
%-----------------------------------------
% input data for nodal coordinate values
%-----------------------------------------
elemsize=1.0/nel;
for i=1:nnode
gcoord(i)=elemsize*(i-1);
end
%-----------------------------------------------------
% input data for nodal connectivity for each element
%-----------------------------------------------------
for i=1:nel
nodes(i,1)=i;
nodes(i,2)=i+1;
end
%-------------------------------------
% input data for boundary conditions
%-------------------------------------
bcdof(1)=1; % first node is constrained
bcval(1)=0; % whose described value is 0
bcdof(2)=nnode; % 4th node is constrained
bcval(2)=1; % whose described value is 0
%-----------------------------------------
% loop for iteration
%-----------------------------------------
error=1; % error is set to 1 arbitrarily
solold=gcoord'; % assume a linear function initially
it=0;
while error > toler
it=it+1; % iteration counter
%-----------------------------------------
% initialization of matrices and vectors
%-----------------------------------------
ff=zeros(sdof,1); % initialization of system force vector
kk=zeros(sdof,sdof); % initialization of system matrix
index=zeros(nnel*ndof,1); % initialization of index vector
%-----------------------------------------------------------------
% computation of element matrices and vectors and their assembly
%-----------------------------------------------------------------
for iel=1:nel % loop for the total number of elements
nl=nodes(iel,1); nr=nodes(iel,2); % extract nodes for (iel)-th element
xl=gcoord(nl); xr=gcoord(nr);% extract nodal coord values for the element
eleng=xr-xl; % element length
index=feeldof1(iel,nnel,ndof);% extract system dofs associated with element
soll=solold(nl); solr=solold(nr); % extract old solutions
%---------------------------------------------------
% element matrix and vector for quasilinearization
%---------------------------------------------------
k(1,1)=1/eleng-alpa*(2*soll+solr)/6+alpa*(solr-soll)/3; % element matrix
k(1,2)=-1/eleng+alpa*(2*soll+solr)/6+alpa*(solr-soll)/6;
k(2,1)=-1/eleng-alpa*(soll+2*solr)/6+alpa*(solr-soll)/6;
k(2,2)=1/eleng+alpa*(soll+2*solr)/6+alpa*(solr-soll)/3;
f(1)=alpa*(solr-soll)*(2*soll+solr)/6; % element vector
f(2)=alpa*(solr-soll)*(soll+2*solr)/6;
[kk,ff]=feasmbl2(kk,ff,k,f,index); % assemble element matrices and vectors
end % end of loop for elments
%-----------------------------
% apply boundary conditions
%-----------------------------
[kk,ff]=feaplyc2(kk,ff,bcdof,bcval);
%----------------------------
% solve the matrix equation
%----------------------------
fsol=kk\ff;
%----------------------------
% check the error
%----------------------------
error=0;
for i=2:(nnode-1)
error=error+(fsol(i)-solold(i))^2/fsol(i)^2;
end
error=sqrt(error);
solold=fsol; % assigne previous solution
end % end of iteration loop
%--------------------------------
% plot of the solution
%--------------------------------
No_of_Iteration=it % print number of iteration for convergency
plot(gcoord',fsol);
xlabel('x-axis')
ylabel('Solution')
%------------------------------