-
Notifications
You must be signed in to change notification settings - Fork 3
/
solver_u.m
506 lines (459 loc) · 32 KB
/
solver_u.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
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
function [u_s] = solver_u(visc_s, visc, AGlen_s)
global rho g n epsilon type_BBC beta2_s lambda_max m_max epr isFlowband...
xi Ms N dx dzeta H_s W W_s dzetadx dzetadx_s dhSdx_s dhBdx_s iter_u...
u_s_lst At_CTS iTimeStep
L1 = 4*visc_s;
L2 = 4*visc_s.*dzetadx_s.^2 + visc_s./(ones(N,1)*H_s).^2;
L3 = 8*dzetadx_s.*visc_s;
L4 = zeros(N,Ms);
L5 = zeros(N,Ms);
L6 = zeros(N,Ms);
Ls1 = zeros(1,Ms);
Ls2 = zeros(1,Ms);
LsW = zeros(1,Ms);
Lb_CFL = zeros(1,Ms);
frozen = zeros(1,Ms);
LHS = zeros(N*Ms, N*Ms);
RHS = zeros(N*Ms,1);
for i = 2:Ms-1
for j = 2:N-1
k = (i-1)*N + j;
RHS(k,1) = rho*g*dhSdx_s(i);
if i == 2
L4(j,i) = 4*visc_s(j,i)*(-4*dzetadx_s(j,i-1)+3*dzetadx_s(j,i)+dzetadx_s(j,i+1))/(3*dx)+ ...
4*visc_s(j,i)*dzetadx_s(j,i)*(dzetadx_s(j+1,i)-dzetadx_s(j-1,i))/(2*dzeta) + ...
4*dzetadx_s(j,i)*(-4*visc_s(j,i-1)+3*visc_s(j,i)+visc_s(j,i+1))/(3*dx) + ...
(visc_s(j+1,i)-visc_s(j-1,i))/(2*dzeta)*(4*dzetadx_s(j,i)^2+1/H_s(i)^2) + ...
dzetadx_s(j,i)*(2*visc_s(j,i)/W_s(j,i))*...
((-4*W_s(j,i-1)+3*W_s(j,i)+W_s(j,i+1))/(3*dx)+dzetadx_s(j,i)*(W_s(j+1,i)-W_s(j-1,i))/(2*dzeta));
L4(N,i) = 4*visc_s(N,i)*(-4*dzetadx_s(N,i-1)+3*dzetadx_s(N,i)+dzetadx_s(N,i+1))/(3*dx)+ ...
4*visc_s(N,i)*dzetadx_s(N,i)*(3*dzetadx_s(N,i)-4*dzetadx_s(N-1,i)+dzetadx_s(N-2,i))/(2*dzeta) + ...
4*dzetadx_s(N,i)*(-4*visc_s(N,i-1)+3*visc_s(N,i)+visc_s(N,i+1))/(3*dx) + ...
(3*visc_s(N,i)-4*visc_s(N-1,i)+visc_s(N-2,i))/(2*dzeta)*(4*dzetadx_s(N,i)^2+1/H_s(i)^2) + ...
dzetadx_s(N,i)*(2*visc_s(N,i)/W_s(N,i))*...
((-4*W_s(N,i-1)+3*W_s(N,i)+W_s(N,i+1))/(3*dx)+dzetadx_s(N,i)*(3*W_s(N,i)-4*W_s(N-1,i)+W_s(N-2,i))/(2*dzeta));
L4(1,i) = 4*visc_s(1,i)*(-4*dzetadx_s(1,i-1)+3*dzetadx_s(1,i)+dzetadx_s(1,i+1))/(3*dx) + ...
4*visc_s(1,i)*dzetadx_s(1,i)*(-3*dzetadx_s(1,i)+4*dzetadx_s(2,i)-dzetadx_s(3,i))/(2*dzeta) + ...
4*dzetadx_s(1,i)*(-4*visc_s(1,i-1)+3*visc_s(1,i)+visc_s(1,i+1))/(3*dx) + ...
(-3*visc_s(1,i)+4*visc_s(2,i)-visc_s(3,i))/(2*dzeta)*(4*dzetadx_s(1,i)^2+1/H_s(i)^2) + ...
dzetadx_s(1,i)*(2*visc_s(1,i)/W_s(1,i))*...
((-4*W_s(1,i-1)+3*W_s(1,i)+W_s(1,i+1))/(3*dx)+dzetadx_s(1,i)*(-3*W_s(1,i)+4*W_s(2,i)-W_s(3,i))/(2*dzeta));
elseif i == Ms-1
L4(j,i) = 4*visc_s(j,i)*(-1*dzetadx_s(j,i-1)-3*dzetadx_s(j,i)+4*dzetadx_s(j,i+1))/(3*dx)+...
4*visc_s(j,i)*dzetadx_s(j,i)*(dzetadx_s(j+1,i)-dzetadx_s(j-1,i))/(2*dzeta) +...
4*dzetadx_s(j,i)*(-1*visc_s(j,i-1)-3*visc_s(j,i)+4*visc_s(j,i+1))/(3*dx) +...
(visc_s(j+1,i)-visc_s(j-1,i))/(2*dzeta)*(4*dzetadx_s(j,i)^2+1/H_s(i)^2) + ...
dzetadx_s(j,i)*(2*visc_s(j,i)/W_s(j,i))*...
((-1*W_s(j,i-1)-3*W_s(j,i)+4*W_s(j,i+1))/(3*dx)+dzetadx_s(j,i)*(W_s(j+1,i)-W_s(j-1,i))/(2*dzeta));
L4(N,i) = 4*visc_s(N,i)*(-1*dzetadx_s(N,i-1)-3*dzetadx_s(N,i)+4*dzetadx_s(N,i+1))/(3*dx)+...
4*visc_s(N,i)*dzetadx_s(N,i)*(3*dzetadx_s(N,i)-4*dzetadx_s(N-1,i)+dzetadx_s(N-2,i))/(2*dzeta) +...
4*dzetadx_s(N,i)*(-1*visc_s(N,i-1)-3*visc_s(N,i)+4*visc_s(N,i+1))/(3*dx) +...
(3*visc_s(N,i)-4*visc_s(N-1,i)+visc_s(N-2,i))/(2*dzeta)*(4*dzetadx_s(N,i)^2+1/H_s(i)^2) + ...
dzetadx_s(N,i)*(2*visc_s(N,i)/W_s(N,i))*...
((-1*W_s(N,i-1)-3*W_s(N,i)+4*W_s(N,i+1))/(3*dx)+dzetadx_s(N,i)*(3*W_s(N,i)-4*W_s(N-1,i)+W_s(N-2,i))/(2*dzeta));
L4(1,i) = 4*visc_s(1,i)*(-1*dzetadx_s(1,i-1)-3*dzetadx_s(1,i)+4*dzetadx_s(1,i+1))/(3*dx)+...
4*visc_s(1,i)*dzetadx_s(1,i)*(-3*dzetadx_s(1,i)+4*dzetadx_s(2,i)-dzetadx_s(3,i))/(2*dzeta) +...
4*dzetadx_s(1,i)*(-1*visc_s(1,i-1)-3*visc_s(1,i)+4*visc_s(1,i+1))/(3*dx) +...
(-3*visc_s(1,i)+4*visc_s(2,i)-visc_s(3,i))/(2*dzeta)*(4*dzetadx_s(1,i)^2+1/H_s(i)^2) + ...
dzetadx_s(1,i)*(2*visc_s(1,i)/W_s(1,i))*...
((-1*W_s(1,i-1)-3*W_s(1,i)+4*W_s(1,i+1))/(3*dx)+dzetadx_s(1,i)*(-3*W_s(1,i)+4*W_s(2,i)-W_s(3,i))/(2*dzeta));
else
L4(j,i) = 4*visc_s(j,i)*(dzetadx(j,i)-dzetadx(j,i-1))/(dx)+...
4*visc_s(j,i)*dzetadx_s(j,i)*(dzetadx_s(j+1,i)-dzetadx_s(j-1,i))/(2*dzeta) +...
4*dzetadx_s(j,i)*(visc(j,i)-visc(j,i-1))/(dx) +...
(visc_s(j+1,i)-visc_s(j-1,i))/(2*dzeta)*(4*dzetadx_s(j,i)^2+1/H_s(i)^2) + ...
dzetadx_s(j,i)*(2*visc_s(j,i)/W_s(j,i))*...
((W(j,i)-W(j,i-1))/(dx)+dzetadx_s(j,i)*(W_s(j+1,i)-W_s(j-1,i))/(2*dzeta));
L4(N,i) = 4*visc_s(N,i)*(dzetadx(N,i)-dzetadx(N,i-1))/(dx)+...
4*visc_s(N,i)*dzetadx_s(N,i)*(3*dzetadx_s(N,i)-4*dzetadx_s(N-1,i)+dzetadx_s(N-2,i))/(2*dzeta) +...
4*dzetadx_s(N,i)*(visc(N,i)-visc(N,i-1))/(dx) +...
(3*visc_s(N,i)-4*visc_s(N-1,i)+visc_s(N-2,i))/(2*dzeta)*(4*dzetadx_s(N,i)^2+1/H_s(i)^2) + ...
dzetadx_s(N,i)*(2*visc_s(N,i)/W_s(N,i))*...
((W(N,i)-W(N,i-1))/(dx)+dzetadx_s(N,i)*(3*W_s(N,i)-4*W_s(N-1,i)+W_s(N-2,i))/(2*dzeta));
L4(1,i) = 4*visc_s(1,i)*(dzetadx(1,i)-dzetadx(1,i-1))/(dx)+...
4*visc_s(1,i)*dzetadx_s(1,i)*(-3*dzetadx_s(1,i)+4*dzetadx_s(2,i)-dzetadx_s(3,i))/(2*dzeta) +...
4*dzetadx_s(1,i)*(visc(1,i)-visc(1,i-1))/(dx) +...
(-3*visc_s(1,i)+4*visc_s(2,i)-visc_s(3,i))/(2*dzeta)*(4*dzetadx_s(1,i)^2+1/H_s(i)^2) + ...
dzetadx_s(1,i)*(2*visc_s(1,i)/W_s(1,i))*...
((W(1,i)-W(1,i-1))/(dx)+dzetadx_s(1,i)*(-3*W_s(1,i)+4*W_s(2,i)-W_s(3,i))/(2*dzeta));
end
if i == 2
L5(j,i) = 4*(-4*visc_s(j,i-1)+3*visc_s(j,i)+visc_s(j,i+1))/(3*dx) +...
4*dzetadx_s(j,i)*(visc_s(j+1,i)-visc_s(j-1,i))/(2*dzeta) + ...
2*visc_s(j,i)/W_s(j,i)*...
((-4*W_s(j,i-1)+3*W_s(j,i)+W_s(j,i+1))/(3*dx)+dzetadx_s(j,i)*(W_s(j+1,i)-W_s(j-1,i))/(2*dzeta));
L5(N,i) = 4*(-4*visc_s(N,i-1)+3*visc_s(N,i)+visc_s(N,i+1))/(3*dx) +...
4*dzetadx_s(N,i)*(3*visc_s(N,i)-4*visc_s(N-1,i)+visc_s(N-2,i))/(2*dzeta) + ...
2*visc_s(N,i)/W_s(N,i)*...
((-4*W_s(N,i-1)+3*W_s(N,i)+W_s(N,i+1))/(3*dx)+dzetadx_s(N,i)*(3*W_s(N,i)-4*W_s(N-1,i)+W_s(N-2,i))/(2*dzeta));
L5(1,i) = 4*(-4*visc_s(1,i-1)+3*visc_s(1,i)+visc_s(1,i+1))/(3*dx) +...
4*dzetadx_s(1,i)*(-3*visc_s(1,i)+4*visc_s(2,i)-visc_s(3,i))/(2*dzeta) + ...
2*visc_s(1,i)/W_s(1,i)*...
((-4*W_s(1,i-1)+3*W_s(1,i)+W_s(1,i+1))/(3*dx)+dzetadx_s(1,i)*(-3*W_s(1,i)+4*W_s(2,i)-W_s(3,i))/(2*dzeta));
elseif i == Ms-1
L5(j,i) = 4*(-1*visc_s(j,i-1)-3*visc_s(j,i)+4*visc_s(j,i+1))/(3*dx) +...
4*dzetadx_s(j,i)*(visc_s(j+1,i)-visc_s(j-1,i))/(2*dzeta) + ...
2*visc_s(j,i)/W_s(j,i)*...
((-1*W_s(j,i-1)-3*W_s(j,i)+4*W_s(j,i+1))/(3*dx)+dzetadx_s(j,i)*(W_s(j+1,i)-W_s(j-1,i))/(2*dzeta));
L5(N,i) = 4*(-1*visc_s(N,i-1)-3*visc_s(N,i)+4*visc_s(N,i+1))/(3*dx) +...
4*dzetadx_s(N,i)*(3*visc_s(N,i)-4*visc_s(N-1,i)+visc_s(N-2,i))/(2*dzeta) + ...
2*visc_s(N,i)/W_s(N,i)*...
((-1*W_s(N,i-1)-3*W_s(N,i)+4*W_s(N,i+1))/(3*dx)+dzetadx_s(N,i)*(3*W_s(N,i)-4*W_s(N-1,i)+W_s(N-2,i))/(2*dzeta));
L5(1,i) = 4*(-1*visc_s(1,i-1)-3*visc_s(1,i)+4*visc_s(1,i+1))/(3*dx) +...
4*dzetadx_s(1,i)*(-3*visc_s(1,i)+4*visc_s(2,i)-visc_s(3,i))/(2*dzeta) + ...
2*visc_s(1,i)/W_s(1,i)*...
((-1*W_s(1,i-1)-3*W_s(1,i)+4*W_s(1,i+1))/(3*dx)+dzetadx_s(1,i)*(-3*W_s(1,i)+4*W_s(2,i)-W_s(3,i))/(2*dzeta));
else
L5(j,i) = 4*(visc(j,i)-visc(j,i-1))/(dx) +...
4*dzetadx_s(j,i)*(visc_s(j+1,i)-visc_s(j-1,i))/(2*dzeta) + ...
2*visc_s(j,i)/W_s(j,i)*...
((W(j,i)-W(j,i-1))/(dx)+dzetadx_s(j,i)*(W_s(j+1,i)-W_s(j-1,i))/(2*dzeta));
L5(N,i) = 4*(visc(N,i)-visc(N,i-1))/(dx) +...
4*dzetadx_s(N,i)*(3*visc_s(N,i)-4*visc_s(N-1,i)+visc_s(N-2,i))/(2*dzeta) + ...
2*visc_s(N,i)/W_s(N,i)*...
((W(N,i)-W(N,i-1))/(dx)+dzetadx_s(N,i)*(3*W_s(N,i)-4*W_s(N-1,i)+W_s(N-2,i))/(2*dzeta));
L5(1,i) = 4*(visc(1,i)-visc(1,i-1))/(dx) +...
4*dzetadx_s(1,i)*(-3*visc_s(1,i)+4*visc_s(2,i)-visc_s(3,i))/(2*dzeta) + ...
2*visc_s(1,i)/W_s(1,i)*...
((W(1,i)-W(1,i-1))/(dx)+dzetadx_s(1,i)*(-3*W_s(1,i)+4*W_s(2,i)-W_s(3,i))/(2*dzeta));
end
if isFlowband
if i == 2
L6(j,i) = 2/W_s(j,i)*((-4*visc_s(j,i-1)+3*visc_s(j,i)+visc_s(j,i+1))/(3*dx)+...
dzetadx_s(j,i)*(visc_s(j+1,i)-visc_s(j-1,i))/(2*dzeta))*...
((-4*W_s(j,i-1)+3*W_s(j,i)+W_s(j,i+1))/(3*dx) + dzetadx_s(j,i)*(W_s(j+1,i)-W_s(j-1,i))/(2*dzeta)) +...
2*visc_s(j,i)/W_s(j,i)*((8*W_s(j,i-1)-12*W_s(j,i)+4*W_s(j,i+1))/(3*dx^2) + ...
dzetadx_s(j,i)^2*(W_s(j-1,i)-2*W_s(j,i)+W_s(j+1,i))/(dzeta^2) + ...
2*dzetadx_s(j,i)*(-4*(W_s(j+1,i-1)-W_s(j-1,i-1))+...
3*(W_s(j+1,i)-W_s(j-1,i))+(W_s(j+1,i+1)-W_s(j-1,i+1)))/(6*dx*dzeta) + ...
(W_s(j+1,i)-W_s(j-1,i))/(2*dzeta)*...
((-4*dzetadx_s(j,i-1)+3*dzetadx_s(j,i)+dzetadx_s(j,i+1))/(3*dx)+...
dzetadx_s(j,i)*(dzetadx_s(j+1,i)-dzetadx_s(j-1,i))/(2*dzeta)) -...
((-4*W_s(j,i-1)+3*W_s(j,i)+W_s(j,i+1))/(3*dx)+dzetadx_s(j,i)*(W_s(j+1,i)-W_s(j-1,i))/(2*dzeta))^2/W_s(j,i)) - ...
visc_s(j,i)/W_s(j,i)^2;
L6(N,i) = 2/W_s(N,i)*((-4*visc_s(N,i-1)+3*visc_s(N,i)+visc_s(N,i+1))/(3*dx)+...
dzetadx_s(N,i)*(3*visc_s(N,i)-4*visc_s(N-1,i)+visc_s(N-2,i))/(2*dzeta))*...
((-4*W_s(N,i-1)+3*W_s(N,i)+W_s(N,i+1))/(3*dx)+dzetadx_s(N,i)*(3*W_s(N,i)-4*W_s(N-1,i)+W_s(N-2,i))/(2*dzeta)) +...
2*visc_s(N,i)/W_s(N,i)*((8*W_s(N,i-1)-12*W_s(N,i)+4*W_s(N,i+1))/(3*dx^2) + ...
dzetadx_s(N,i)^2*(W_s(N,i)-2*W_s(N-1,i)+W_s(N-2,i))/(dzeta^2) + ...
2*dzetadx_s(N,i)*(-4*(3*W_s(N,i-1)-4*W_s(N-1,i-1)+W_s(N-2,i-1))+...
3*(3*W_s(N,i)-4*W_s(N-1,i)+W_s(N-2,i))+(3*W_s(N,i+1)-4*W_s(N-1,i+1)+W_s(N-2,i+1)))/(6*dx*dzeta) + ...
(3*W_s(N,i)-4*W_s(N-1,i)+W_s(N-2,i))/(2*dzeta)*...
((-4*dzetadx_s(N,i-1)+3*dzetadx_s(N,i)+dzetadx_s(N,i+1))/(3*dx) +...
dzetadx_s(N,i)*(3*dzetadx_s(N,i)-4*dzetadx_s(N-1,i)+dzetadx_s(N-2,i))/(2*dzeta)) -...
((-4*W_s(N,i-1)+3*W_s(N,i)+W_s(N,i+1))/(3*dx)+dzetadx_s(N,i)*(3*W_s(N,i)-4*W_s(N-1,i)+W_s(N-2,i))/(2*dzeta))^2/W_s(N,i)) - ...
visc_s(N,i)/W_s(N,i)^2;
L6(1,i) = 2/W_s(1,i)*((-4*visc_s(1,i-1)+3*visc_s(1,i)+visc_s(1,i+1))/(3*dx)+...
dzetadx_s(1,i)*(-3*visc_s(1,i)+4*visc_s(2,i)-visc_s(3,i))/(2*dzeta))*...
((-4*W_s(1,i-1)+3*W_s(1,i)+W_s(1,i+1))/(3*dx)+dzetadx_s(1,i)*(-3*W_s(1,i)+4*W_s(2,i)-W_s(3,i))/(2*dzeta)) +...
2*visc_s(1,i)/W_s(1,i)*((8*W_s(1,i-1)-12*W_s(1,i)+4*W_s(1,i+1))/(3*dx^2) + ...
dzetadx_s(1,i)^2*(W_s(1,i)-2*W_s(2,i)+W_s(3,i))/(dzeta^2) +...
2*dzetadx_s(1,i)*(-4*(-3*W_s(1,i-1)+4*W_s(2,i-1)+W_s(3,i-1))+...
3*(-3*W_s(1,i)+4*W_s(2,i)+W_s(3,i))+(-3*W_s(1,i+1)+4*W_s(2,i+1)+W_s(3,i+1)))/(6*dx*dzeta) + ...
(-3*W_s(1,i)+4*W_s(2,i)-W_s(3,i))/(2*dzeta)*...
((-4*dzetadx_s(1,i-1)+3*dzetadx_s(1,i)+dzetadx_s(1,i+1))/(3*dx)+...
dzetadx_s(1,i)*(-3*dzetadx_s(1,i)+4*dzetadx_s(2,i)-dzetadx_s(3,i))/(2*dzeta)) -...
((-4*W_s(1,i-1)+3*W_s(1,i)+W_s(1,i+1))/(3*dx)+dzetadx_s(1,i)*(-3*W_s(1,i)+4*W_s(2,i)-W_s(3,i))/(2*dzeta))^2/W_s(1,i)) - ...
visc_s(1,i)/W_s(1,i)^2;
elseif i == Ms-1
L6(j,i) = 2/W_s(j,i)*((-1*visc_s(j,i-1)-3*visc_s(j,i)+4*visc_s(j,i+1))/(3*dx)+...
dzetadx_s(j,i)*(visc_s(j+1,i)-visc_s(j-1,i))/(2*dzeta))*...
((-1*W_s(j,i-1)-3*W_s(j,i)+4*W_s(j,i+1))/(3*dx)+dzetadx_s(j,i)*(W_s(j+1,i)-W_s(j-1,i))/(2*dzeta)) +...
2*visc_s(j,i)/W_s(j,i)*((4*W_s(j,i-1)-12*W_s(j,i)+8*W_s(j,i+1))/(3*dx^2) + ...
dzetadx_s(j,i)^2*(W_s(j-1,i)-2*W_s(j,i)+W_s(j+1,i))/(dzeta^2) +...
2*dzetadx_s(j,i)*(-1*(W_s(j-1,i-1)+W_s(j,i-1)+W_s(j+1,i-1))-...
3*(W_s(j-1,i)+W_s(j,i)+W_s(j+1,i))+4*(W_s(j-1,i+1)+W_s(j,i+1)+W_s(j+1,i+1)))/(6*dx*dzeta) +...
(W_s(j+1,i)-W_s(j-1,i))/(2*dzeta)*...
((-1*dzetadx_s(j,i-1)-3*dzetadx_s(j,i)+4*dzetadx_s(j,i+1))/(3*dx)+...
dzetadx_s(j,i)*(dzetadx_s(j+1,i)-dzetadx_s(j-1,i))/(2*dzeta)) -...
((-1*W_s(j,i-1)-3*W_s(j,i)+4*W_s(j,i+1))/(3*dx)+dzetadx_s(j,i)*(W_s(j+1,i)-W_s(j-1,i))/(2*dzeta))^2/W_s(j,i)) - ...
visc_s(j,i)/W_s(j,i)^2;
L6(N,i) = 2/W_s(N,i)*((-1*visc_s(N,i-1)-3*visc_s(N,i)+4*visc_s(N,i+1))/(3*dx)+...
dzetadx_s(N,i)*(3*visc_s(N,i)-4*visc_s(N-1,i)+visc_s(N-2,i))/(2*dzeta))*...
((-1*W_s(N,i-1)-3*W_s(N,i)+4*W_s(N,i+1))/(3*dx)+dzetadx_s(N,i)*(3*W_s(N,i)-4*W_s(N-1,i)+W_s(N-2,i))/(2*dzeta)) +...
2*visc_s(N,i)/W_s(N,i)*((4*W_s(N,i-1)-12*W_s(N,i)+8*W_s(N,i+1))/(3*dx^2) + ...
dzetadx_s(N,i)^2*(W_s(N,i)-2*W_s(N-1,i)+W_s(N-2,i))/(dzeta^2) +...
2*dzetadx_s(N,i)*(-1*(3*W_s(N,i-1)-4*W_s(N-1,i-1)+W_s(N,i-1))-...
3*(3*W_s(N,i)-4*W_s(N-1,i)+W_s(N,i))+4*(3*W_s(N,i+1)-4*W_s(N-1,i+1)+W_s(N,i+1)))/(6*dx*dzeta) + ...
(3*W_s(N,i)-4*W_s(N-1,i)+W_s(N-2,i))/(2*dzeta)*...
((-1*dzetadx_s(N,i-1)-3*dzetadx_s(N,i)+4*dzetadx_s(N,i+1))/(3*dx)+...
dzetadx_s(N,i)*(3*dzetadx_s(N,i)-4*dzetadx_s(N-1,i)+dzetadx_s(N-2,i))/(2*dzeta)) -...
((-1*W_s(N,i-1)-3*W_s(N,i)+4*W_s(N,i+1))/(3*dx)+dzetadx_s(N,i)*(3*W_s(N,i)-4*W_s(N-1,i)+W_s(N-2,i))/(2*dzeta))^2/W_s(N,i)) - ...
visc_s(N,i)/W_s(N,i)^2;
L6(1,i) = 2/W_s(1,i)*((-1*visc_s(1,i-1)-3*visc_s(1,i)+4*visc_s(1,i+1))/(3*dx)+...
dzetadx_s(1,i)*(-3*visc_s(1,i)+4*visc_s(2,i)-visc_s(3,i))/(2*dzeta))*...
((-1*W_s(1,i-1)-3*W_s(1,i)+4*W_s(1,i+1))/(3*dx)+dzetadx_s(1,i)*(-3*W_s(1,i)+4*W_s(2,i)-W_s(3,i))/(2*dzeta)) + ...
2*visc_s(1,i)/W_s(1,i)*((4*W_s(1,i-1)-12*W_s(1,i)+8*W_s(1,i+1))/(3*dx^2) + ...
dzetadx_s(1,i)^2*(W_s(1,i)-2*W_s(2,i)+W_s(3,i))/(dzeta^2) + ...
2*dzetadx_s(1,i)*(-1*(-3*W_s(1,i-1)+4*W_s(2,i-1)+W_s(3,i-1))-...
3*(-3*W_s(1,i)+4*W_s(2,i)+W_s(3,i))+4*(-3*W_s(1,i+1)+4*W_s(2,i+1)+W_s(3,i+1)))/(6*dx*dzeta) + ...
(-3*W_s(1,i)+4*W_s(2,i)-W_s(3,i))/(2*dzeta)*...
((-1*dzetadx_s(1,i-1)-3*dzetadx_s(1,i)+4*dzetadx_s(1,i+1))/(3*dx)+...
dzetadx_s(1,i)*(-3*dzetadx_s(1,i)+4*dzetadx_s(2,i)-dzetadx_s(3,i))/(2*dzeta)) - ...
((-1*W_s(1,i-1)-3*W_s(1,i)+4*W_s(1,i+1))/(3*dx)+dzetadx_s(1,i)*(-3*W_s(1,i)+4*W_s(2,i)-W_s(3,i))/(2*dzeta))^2/W_s(1,i)) - ...
visc_s(1,i)/W_s(1,i)^2;
else
L6(j,i) = 2/W_s(j,i)*((visc(j,i)-visc(j,i-1))/(dx) + dzetadx_s(j,i)*(visc_s(j+1,i)-visc_s(j-1,i))/(2*dzeta))*...
((W(j,i)-W(j,i-1))/(dx) + dzetadx_s(j,i)*(W_s(j+1,i)-W_s(j-1,i))/(2*dzeta)) +...
2*visc_s(j,i)/W_s(j,i)*((W_s(j,i-1)-2*W_s(j,i)+W_s(j,i+1))/(dx^2) + ...
dzetadx_s(j,i)^2*(W_s(j-1,i)-2*W_s(j,i)+W_s(j+1,i))/(dzeta^2) +...
2*dzetadx_s(j,i)*(W_s(j+1,i+1)+W_s(j-1,i-1)-W_s(j+1,i-1)-W_s(j-1,i+1))/(4*dx*dzeta) + ...
(W_s(j+1,i)-W_s(j-1,i))/(2*dzeta)*...
((dzetadx(j,i)-dzetadx(j,i-1))/(dx)+dzetadx_s(j,i)*(dzetadx_s(j+1,i)-dzetadx_s(j-1,i))/(2*dzeta)) -...
((W(j,i)-W(j,i-1))/(dx)+dzetadx_s(j,i)*(W_s(j+1,i)-W_s(j-1,i))/(2*dzeta))^2/W_s(j,i)) - ...
visc_s(j,i)/W_s(j,i)^2;
L6(N,i) = 2/W_s(N,i)*((visc(N,i)-visc(N,i-1))/(dx) + dzetadx_s(N,i)*(3*visc_s(N,i)-4*visc_s(N-1,i)+visc_s(N-2,i))/(2*dzeta))*...
((W(N,i)-W(N,i-1))/(dx)+dzetadx_s(N,i)*(3*W_s(N,i)-4*W_s(N-1,i)+W_s(N-2,i))/(2*dzeta)) +...
2*visc_s(N,i)/W_s(N,i)*((W_s(N,i-1)-2*W_s(N,i)+W_s(N,i+1))/(dx^2) + ...
dzetadx_s(N,i)^2*(W_s(N,i)-2*W_s(N-1,i)+W_s(N-2,i))/(dzeta^2) + ...
2*dzetadx_s(N,i)*((3*W_s(N,i+1)-4*W_s(N-1,i+1)+W_s(N-2,i+1))-...
(3*W_s(N,i-1)-4*W_s(N-1,i-1)+W_s(N-2,i-1)))/(4*dx*dzeta) + ...
(3*W_s(N,i)-4*W_s(N-1,i)+W_s(N-2,i))/(2*dzeta)*...
((dzetadx(N,i)-dzetadx(N,i-1))/(dx)+dzetadx_s(N,i)*(3*dzetadx_s(N,i)-4*dzetadx_s(N-1,i)+dzetadx_s(N-2,i))/(2*dzeta)) -...
((W(N,i)-W(N,i-1))/(dx)+dzetadx_s(N,i)*(3*W_s(N,i)-4*W_s(N-1,i)+W_s(N-2,i))/(2*dzeta))^2/W_s(N,i)) - ...
visc_s(N,i)/W_s(N,i)^2;
L6(1,i) = 2/W_s(1,i)*((visc(1,i)-visc(1,i-1))/(dx)+dzetadx_s(1,i)*(-3*visc_s(1,i)+4*visc_s(2,i)-visc_s(3,i))/(2*dzeta))*...
((W(1,i)-W(1,i-1))/(dx)+dzetadx_s(1,i)*(-3*W_s(1,i)+4*W_s(2,i)-W_s(3,i))/(2*dzeta)) +...
2*visc_s(1,i)/W_s(1,i)*((W_s(1,i-1)-2*W_s(1,i)+W_s(1,i+1))/(dx^2) + dzetadx_s(1,i)^2*(W_s(1,i)-2*W_s(2,i)+W_s(3,i))/(dzeta^2) +...
2*dzetadx_s(1,i)*((-3*W_s(1,i+1)+4*W_s(2,i+1)-W_s(3,i+1))-(-3*W_s(1,i-1)+4*W_s(2,i-1)-W_s(3,i-1)))/(2*dx*dzeta) + ...
(-3*W_s(1,i)+4*W_s(2,i)-W_s(3,i))/(2*dzeta)*...
((dzetadx(1,i)-dzetadx(1,i-1))/(dx)+dzetadx_s(1,i)*(-3*dzetadx_s(1,i)+4*dzetadx_s(2,i)-dzetadx_s(3,i))/(2*dzeta)) -...
((W(1,i)-W(1,i-1))/(dx)+dzetadx_s(1,i)*(-3*W_s(1,i)+4*W_s(2,i)-W_s(3,i))/(2*dzeta))^2/W_s(1,i)) - visc_s(1,i)/W_s(1,i)^2;
end
end
if i == 2
LHS(k,k-N-1) = 2*L3(j,i)/(3*dx*dzeta);
LHS(k,k-N) = 8*L1(j,i)/(3*dx^2)-4*L5(j,i)/(3*dx);
LHS(k,k-N+1) = -2*L3(j,i)/(3*dx*dzeta);
LHS(k,k-1) = L2(j,i)/dzeta^2-L3(j,i)/(2*dx*dzeta)-L4(j,i)/(2*dzeta);
LHS(k,k) = -4*L1(j,i)/(dx^2)-2*L2(j,i)/dzeta^2+L5(j,i)/dx+L6(j,i);
LHS(k,k+1) = L2(j,i)/(dzeta^2)+L3(j,i)/(2*dx*dzeta)+L4(j,i)/(2*dzeta);
LHS(k,k+N-1) = -L3(j,i)/(6*dx*dzeta);
LHS(k,k+N) = 4*L1(j,i)/(3*dx^2)+L5(j,i)/(3*dx);
LHS(k,k+N+1) = L3(j,i)/(6*dx*dzeta);
elseif i == Ms-1
LHS(k,k-N-1) = L3(j,i)/(6*dx*dzeta);
LHS(k,k-N) = 4*L1(j,i)/(3*dx^2)-L5(j,i)/(3*dx);
LHS(k,k-N+1) = -L3(j,i)/(6*dx*dzeta);
LHS(k,k-1) = L2(j,i)/(dzeta^2)+L3(j,i)/(2*dx*dzeta)-L4(j,i)/(2*dzeta);
LHS(k,k) = -4*L1(j,i)/(dx^2)-2*L2(j,i)/dzeta^2-L5(j,i)/dx+L6(j,i);
LHS(k,k+1) = L2(j,i)/(dzeta^2)-L3(j,i)/(2*dx*dzeta)+L4(j,i)/(2*dzeta);
LHS(k,k+N-1) = -2*L3(j,i)/(3*dx*dzeta);
LHS(k,k+N) = 8*L1(j,i)/(3*dx^2)+4*L5(j,i)/(3*dx);
LHS(k,k+N+1) = 2*L3(j,i)/(3*dx*dzeta);
else
LHS(k,k-N-1) = L3(j,i)/(4*dx*dzeta);
LHS(k,k-N) = L1(j,i)/(dx^2) - L5(j,i)/(2*dx);
LHS(k,k-N+1) = -L3(j,i)/(4*dx*dzeta);
LHS(k,k-1) = L2(j,i)/(dzeta^2) - L4(j,i)/(2*dzeta);
LHS(k,k) = -2*L1(j,i)/(dx^2) - 2*L2(j,i)/dzeta^2 + L6(j,i);
LHS(k,k+1) = L2(j,i)/dzeta^2 + L4(j,i)/(2*dzeta);
LHS(k,k+N-1) = -L3(j,i)/(4*dx*dzeta);
LHS(k,k+N) = L1(j,i)/(dx^2) + L5(j,i)/(2*dx);
LHS(k,k+N+1) = L3(j,i)/(4*dx*dzeta);
end
end
end
for i = 1:Ms
RHS((i-1)*N+N,1) = rho*g*dhSdx_s(i);
Ls1(i) = 4*H_s(i)*dhSdx_s(i)/(1-4*H_s(i)*dzetadx_s(N,i)*dhSdx_s(i))*dzeta/dx;
Ls2(i) = 4*H_s(i)*dhSdx_s(i)/(1-4*H_s(i)*dzetadx_s(N,i)*dhSdx_s(i))*dzeta/W_s(i);
if i == 1
LsW(i) = (-8*W_s(N,i)+9*W_s(N,i+1)-W_s(N,i+2))/(3*dx) + dzetadx_s(N,i)*(3*W_s(N,i)-4*W_s(N-1,i)+W_s(N-2,i))/(2*dzeta);
elseif i == 2
LsW(i) = (-4*W_s(N,i-1)+3*W_s(N,i)+W_s(N,i+1))/(3*dx) + dzetadx_s(N,i)*(3*W_s(N,i)-4*W_s(N-1,i)+W_s(N-2,i))/(2*dzeta);
elseif i == Ms-1
LsW(i) = (-W_s(N,i-1)-3*W_s(N,i)+4*W_s(N,i+1))/(3*dx) + dzetadx_s(N,i)*(3*W_s(N,i)-4*W_s(N-1,i)+W_s(N-2,i))/(2*dzeta);
elseif i == Ms
LsW(i) = (8*W_s(N,i)-9*W_s(N,i-1)+W_s(N,i-2))/(3*dx) + dzetadx_s(N,i)*(3*W_s(N,i)-4*W_s(N-1,i)+W_s(N-2,i))/(2*dzeta);
else
LsW(i) = (W_s(N,i+1)-W_s(N,i-1))/(2*dx) + dzetadx_s(N,i)*(3*W_s(N,i)-4*W_s(N-1,i)+W_s(N-2,i))/(2*dzeta);
end
end
LsW = epsilon*LsW;
for i = 2:Ms-1
if i == 2
LHS(i*N, (i-1)*N) = 8*L1(N,i)/(3*dx^2) - 8/3*L2(N,i)*Ls1(i)/(dzeta^2) + ...
L3(N,i)*(64/3*Ls1(i-1)-4*Ls2(i-1)*LsW(i-1)-8*Ls1(i))/(6*dx*dzeta) - ...
8/3*L4(N,i)*Ls1(i)/(2*dzeta) - 4*L5(N,i)/(3*dx);
LHS(i*N, i*N-1) = 2*L2(N,i)/(dzeta^2);
LHS(i*N, i*N) = -12*L1(N,i)/(3*dx^2) + L2(N,i)*(2*Ls1(i)+Ls2(i)*LsW(i)-2)/(dzeta^2) + ...
L3(N,i)*(-24*Ls1(i-1)+6*Ls1(i)+Ls2(i)*LsW(i)-Ls1(i+1))/(6*dx*dzeta) + ...
L4(N,i)*(2*Ls1(i)+Ls2(i)*LsW(i))/(2*dzeta) + 3*L5(N,i)/(3*dx);
LHS(i*N, (i+1)*N) = 4*L1(N,i)/(3*dx^2) + 2/3*L2(N,i)*Ls1(i)/(dzeta^2) + ...
L3(N,i)*(8/3*Ls1(i-1)+2*Ls1(i)+Ls2(i+1)*LsW(i+1))/(6*dx*dzeta) + ...
L4(N,i)*(2/3*Ls1(i))/(2*dzeta) + L5(N,i)/(3*dx);
LHS(i*N, (i+2)*N) = L3(N,i)*Ls1(i+1)/(6*dx*dzeta);
elseif i == Ms-1
LHS(i*N,(i-2)*N) = L3(N,i)*Ls1(i-1)/(6*dx*dzeta);
LHS(i*N,(i-1)*N) = 4*L1(N,i)/(3*dx^2) - 2/3*L2(N,i)*Ls1(i)/(dzeta^2) + ...
L3(N,i)*(-Ls2(i-1)*LsW(i-1)+2*Ls1(i)+8/3*Ls1(i+1))/(6*dx*dzeta)...
-2/3*L4(N,i)*Ls1(i)/(2*dzeta) - L5(N,i)/(3*dx);
LHS(i*N,i*N-1) = 2*L2(N,i)/(dzeta^2);
LHS(i*N,i*N) = -12*L1(N,i)/(3*dx^2) + L2(N,i)*(-2*Ls1(i)+Ls2(i)*LsW(i)-2)/(dzeta^2) + ...
L3(N,i)*(-Ls1(i-1)+6*Ls1(i)-3*Ls2(i)*LsW(i)-24*Ls1(i+1))/(6*dx*dzeta) + ...
L4(N,i)*(-2*Ls1(i)+Ls2(i)*LsW(i))/(2*dzeta) - 3*L5(N,i)/(3*dx);
LHS(i*N,(i+1)*N) = 8*L1(N,i)/(3*dx^2) + 8/3*L2(N,i)*Ls1(i)/(dzeta^2) + ...
L3(N,i)*(-8*Ls1(i)+64/3*Ls1(i+1)+4*Ls2(i+1)*LsW(i+1))/(6*dx*dzeta) + ...
8/3*L4(N,i)*Ls1(i)/(2*dzeta) + 4*L5(N,i)/(3*dx);
else
LHS(i*N,(i-2)*N) = L3(N,i)*Ls1(i-1)/(4*dx*dzeta);
LHS(i*N,(i-1)*N) = L1(N,i)/(dx^2) - L2(N,i)*Ls1(i)/(dzeta^2) - ...
L3(N,i)*Ls2(i-1)*LsW(i-1)/(4*dx*dzeta) - L4(N,i)*Ls1(i)/(2*dzeta) - L5(N,i)/(2*dx);
LHS(i*N,i*N-1) = 2*L2(N,i)/(dzeta^2);
LHS(i*N,i*N) = -2*L1(N,i)/(dx^2) + L2(N,i)*(Ls2(i)*LsW(i)-2)/(dzeta^2) -...
L3(N,i)*(Ls1(i-1)+Ls1(i+1))/(4*dx*dzeta) + L4(N,i)*Ls2(i)*LsW(i)/(2*dzeta);
LHS(i*N,(i+1)*N) = L1(N,i)/(dx^2) + L2(N,i)*Ls1(i)/(dzeta^2) + L3(N,i)*Ls2(i+1)*LsW(i+1)/(4*dx*dzeta) +...
L4(N,i)*Ls1(i)/(2*dzeta) + L5(N,i)/(2*dx);
LHS(i*N,(i+2)*N) = L3(N,i)*Ls1(i+1)/(4*dx*dzeta);
end
end
switch type_BBC
case 1
for i = 2:Ms-1
LHS((i-1)*N+1,(i-1)*N+1) = 1;
RHS((i-1)*N+1,1) = 0;
end
case 2
Cb = 0.84*m_max;
Neff = epr.*(rho*g*H_s);
for i=1:Ms
if iter_u==1
ubi = 0;
else
ubi = u_s_lst(1,i);
end
% Lb_CFL(i) = 4*AGlen_s(1,i)*H_s(i)*dzeta*Cb(i)^n*Neff(i)^n*m_max(i)/...
% (ubi*m_max(i) + Cb(i)^n*Neff(i)^n*lambda_max(i)*AGlen_s(1,i)); % /(1+(H_s(i)/W_s(1,i))^2)
Lb_CFL(i) = 4*AGlen_s(1,i)*H_s(i)*dzeta*Cb^n*Neff(i)^n*m_max/...
(ubi*m_max + Cb^n*Neff(i)^n*lambda_max*AGlen_s(1,i)); % /(1+(H_s(i)/W_s(1,i))^2)
end
frozeni = 0;
for i=2:Ms-1
if iTimeStep==1
frozeni = frozeni + 1;
frozen(frozeni) = i;
LHS((i-1)*N+1,(i-1)*N+1) = 1;
RHS((i-1)*N+1,1) = 0;
else
if At_CTS(iTimeStep-1,i)~=0
if i == 2
LHS((i-1)*N+1,(i-2)*N+1) = 8*L1(1,i)/(3*dx^2) - 4*Lb_CFL(i-1)*L3(1,i)/(6*dx*dzeta) - 4*L5(1,i)/(3*dx);
LHS((i-1)*N+1,(i-1)*N+1) = -12*L1(1,i)/(3*dx^2) - (2+Lb_CFL(i))*L2(1,i)/(dzeta^2) +...
3*L3(1,i)*Lb_CFL(i)/(6*dx*dzeta) + Lb_CFL(i)*L4(1,i)/(2*dzeta) + 3*L5(1,i)/(3*dx) + L6(1,i);
LHS((i-1)*N+1,(i-1)*N+2) = 2*L2(1,i)/dzeta^2;
LHS((i-1)*N+1,i*N+1) = 4*L1(1,i)/(3*dx^2) + Lb_CFL(i+1)*L3(1,i)/(6*dx*dzeta) + L5(1,i)/(3*dx);
elseif i == Ms-1
LHS((i-1)*N+1,(i-2)*N+1) = 4*L1(1,i)/(3*dx^2) - Lb_CFL(i-1)*L3(1,i)/(6*dx*dzeta) - L5(1,i)/(3*dx);
LHS((i-1)*N+1,(i-1)*N+1) = -12*L1(1,i)/(3*dx^2) - (2+Lb_CFL(i))*L2(1,i)/(dzeta^2) -...
3*Lb_CFL(i)*L3(1,i)/(6*dx*dzeta) + Lb_CFL(i)*L4(1,i)/(2*dzeta) - 3*L5(1,i)/(3*dx) + L6(1,i);
LHS((i-1)*N+1,(i-1)*N+2) = 2*L2(1,i)/dzeta^2;
LHS((i-1)*N+1,i*N+1) = 8*L1(1,i)/(3*dx^2) + 4*Lb_CFL(i+1)*L3(1,i)/(6*dx*dzeta) + 4*L5(1,i)/(3*dx);
else
LHS((i-1)*N+1,(i-2)*N+1) = L1(1,i)/(dx^2) - Lb_CFL(i-1)*L3(1,i)/(4*dx*dzeta) - L5(1,i)/(2*dx);
LHS((i-1)*N+1,(i-1)*N+1) = -2*L1(1,i)/(dx^2) - (2+Lb_CFL(i))*L2(1,i)/(dzeta^2) +...
Lb_CFL(i)*L4(1,i)/(2*dzeta) + L6(1,i);
LHS((i-1)*N+1,(i-1)*N+2) = 2*L2(1,i)/(dzeta^2);
LHS((i-1)*N+1,i*N+1) = L1(1,i)/(dx^2) + Lb_CFL(i+1)*L3(1,i)/(4*dx*dzeta) + L5(1,i)/(2*dx);
end
else
frozeni = frozeni + 1;
frozen(frozeni) = i;
LHS((i-1)*N+1,(i-1)*N+1) = 1;
RHS((i-1)*N+1,1) = 0;
end
end
end
case 3
Lb1_LFL = 4.*dzeta.*H_s.*dhBdx_s./(1 - 4.*H_s.*dhBdx_s.*dzetadx_s(1,:))./dx;
Lb2_LFL = 2.*dzeta.*H_s.*beta2_s./visc_s(1,:)./(1 - 4.*H_s.*dhBdx_s.*dzetadx_s(1,:));
for i = 2:Ms-1
RHS((i-1)*N+1,1) = rho*g*dhSdx_s(i);
if i == 2
LHS((i-1)*N+1, (i-2)*N+1) = 8*L1(1,i)/(3*dx^2) + 8*Lb1_LFL(i)*L2(1,i)/(3*dzeta^2) +...
(64*Lb1_LFL(i-1)/3-4*Lb2_LFL(i-1)-8*Lb1_LFL(i))*L3(1,i)/(6*dx*dzeta) -...
8/3*Lb1_LFL(i)*L4(1,i)/(2*dzeta) - 4*L5(1,i)/(3*dx);
LHS((i-1)*N+1, (i-1)*N+1) = -12*L1(1,i)/(3*dx^2) - (2+2*Lb1_LFL(i)+Lb2_LFL(i))*L2(1,i)/(dzeta^2) +...
(-24*Lb1_LFL(i-1)+6*Lb1_LFL(i)+3*Lb2_LFL(i)-Lb1_LFL(i+1))*L3(1,i)/(6*dx*dzeta) +...
(2*Lb1_LFL(i)+Lb2_LFL(i))*L4(1,i)/(2*dzeta) + 3*L5(1,i)/(3*dx) + L6(1,i);
LHS((i-1)*N+1, (i-1)*N+2) = 2*L2(1,i)/(dzeta^2);
LHS((i-1)*N+1, i*N+1) = 4*L1(1,i)/(3*dx^2) - 2/3*Lb1_LFL(i)*L2(1,i)/(dzeta^2) +...
(8/3*Lb1_LFL(i-1)+2*Lb1_LFL(i)+Lb2_LFL(i+1))*L3(1,i)/(6*dx*dzeta) +...
2/3*Lb1_LFL(i)*L4(1,i)/(2*dzeta) + L5(1,i)/(3*dx);
elseif i == Ms-1
LHS((i-1)*N+1, (i-3)*N+1) = Lb1_LFL(i-1)*L3(1,i)/(6*dx*dzeta);
LHS((i-1)*N+1, (i-2)*N+1) = 4*L1(1,i)/(3*dx^2) + 2/3*Lb1_LFL(i)*L2(1,i)/(dzeta^2) +...
(-Lb2_LFL(i-1)+2*Lb1_LFL(i)+8/3*Lb1_LFL(i+1))*L3(1,i)/(6*dx*dzeta) -...
2/3*Lb1_LFL(i)*L4(1,i)/(2*dzeta) - L5(1,i)/(3*dx);
LHS((i-1)*N+1, (i-1)*N+1) = -12*L1(1,i)/(3*dx^2) + (-2+2*Lb1_LFL(i)-Lb2_LFL(i))*L2(1,i)/(dzeta^2) + ...
(-Lb1_LFL(i-1)+6*Lb1_LFL(i)-3*Lb2_LFL(i)-24*Lb1_LFL(i+1))*L3(1,i)/(6*dx*dzeta) +...
(-2*Lb1_LFL(i)+Lb2_LFL(i))*L4(1,i)/(2*dzeta) - 3*L5(1,i)/(3*dx) + L6(1,i);
LHS((i-1)*N+1, (i-1)*N+2) = 2*L2(1,i)/(dzeta^2);
LHS((i-1)*N+1, i*N+1) = 8*L1(1,i)/(3*dx^2) - 8/3*Lb1_LFL(i)*L2(1,i)/(dzeta^2) +...
(-8*Lb1_LFL(i)+64/3*Lb1_LFL(i+1)+4*Lb2_LFL(i+1))*L3(1,i)/(6*dx*dzeta) +...
8/3*Lb1_LFL(i)*L4(1,i)/(2*dzeta) + 4*L5(1,i)/(3*dx);
else
LHS((i-1)*N+1, (i-3)*N+1) = Lb1_LFL(i-1)*L3(1,i)/(4*dx*dzeta);
LHS((i-1)*N+1, (i-2)*N+1) = L1(1,i)/(dx^2) + Lb1_LFL(i)*L2(1,i)/(dzeta^2) - Lb2_LFL(i-1)*L3(1,i)/(4*dx*dzeta) -...
Lb1_LFL(i)*L4(1,i)/(2*dzeta) - L5(1,i)/(2*dx);
LHS((i-1)*N+1, (i-1)*N+1) = -2*L1(1,i)/(dx^2) - (2+Lb2_LFL(i))*L2(1,i)/(dzeta^2) -...
(Lb1_LFL(i+1)+Lb1_LFL(i-1))*L3(1,i)/(4*dx*dzeta) + Lb2_LFL(i)*L4(1,i)/(2*dzeta) + L6(1,i);
LHS((i-1)*N+1, (i-1)*N+2) = 2*L2(1,i)/(dzeta^2);
LHS((i-1)*N+1, i*N+1) = L1(1,i)/(dx^2) - Lb1_LFL(i)*L2(1,i)/(dzeta^2) +...
Lb2_LFL(i+1)*L3(1,i)/(4*dx*dzeta) + Lb1_LFL(i)*L4(1,i)/(2*dzeta) + L5(1,i)/(2*dx);
LHS((i-1)*N+1, (i+1)*N+1) = Lb1_LFL(i+1)*L3(1,i)/(4*dx*dzeta);
end
end
case 4
Lb1_LFL = 4.*dzeta.*H_s.*dhBdx_s./(1 - 4.*H_s.*dhBdx_s.*dzetadx_s(1,:))./dx;
Lb2_LFL = zeros(N,Ms);
xi_s = [xi(1), (xi(1:end-1)+xi(2:end))/2, xi(end)];
for i = 2:Ms-1
if xi_s(i) >= 2200 && xi_s(i) <= 2500
% if xi(i) >= 2200 && xi(i) <= 2500
RHS((i-1)*N+1,1) = rho*g*dhSdx_s(i);
LHS((i-1)*N+1, (i-3)*N+1) = Lb1_LFL(i-1)*L3(1,i)/(4*dx*dzeta);
LHS((i-1)*N+1, (i-2)*N+1) = L1(1,i)/(dx^2) + Lb1_LFL(i)*L2(1,i)/(dzeta^2) - Lb2_LFL(i-1)*L3(1,i)/(4*dx*dzeta) -...
Lb1_LFL(i)*L4(1,i)/(2*dzeta) - L5(1,i)/(2*dx);
LHS((i-1)*N+1, (i-1)*N+1) = -2*L1(1,i)/(dx^2) - (2+Lb2_LFL(i))*L2(1,i)/(dzeta^2) -...
(Lb1_LFL(i+1)+Lb1_LFL(i-1))*L3(1,i)/(4*dx*dzeta) + Lb2_LFL(i)*L4(1,i)/(2*dzeta) + L6(1,i);
LHS((i-1)*N+1, (i-1)*N+2) = 2*L2(1,i)/(dzeta^2);
LHS((i-1)*N+1, i*N+1) = L1(1,i)/(dx^2) - Lb1_LFL(i)*L2(1,i)/(dzeta^2) +...
Lb2_LFL(i+1)*L3(1,i)/(4*dx*dzeta) + Lb1_LFL(i)*L4(1,i)/(2*dzeta) + L5(1,i)/(2*dx);
LHS((i-1)*N+1, (i+1)*N+1) = Lb1_LFL(i+1)*L3(1,i)/(4*dx*dzeta);
else
LHS((i-1)*N+1,(i-1)*N+1) = 1;
RHS((i-1)*N+1,1) = 0;
end
end
case 5
Lb2_LFL = 2.*dzeta.*H_s.*beta2_s./visc_s(1,:);
for i = 2:Ms-1
RHS((i-1)*N+1,1) = rho*g*dhSdx_s(i);
LHS((i-1)*N+1, (i-2)*N+1) = L1(1,i)/(dx^2) - Lb2_LFL(i-1)*L3(1,i)/(4*dx*dzeta) - L5(1,i)/(2*dx);
LHS((i-1)*N+1, (i-1)*N+1) = -2*L1(1,i)/(dx^2) - (2+Lb2_LFL(i))*L2(1,i)/(dzeta^2) +...
Lb2_LFL(i)*L4(1,i)/(2*dzeta) + L6(1,i);
LHS((i-1)*N+1, (i-1)*N+2) = 2*L2(1,i)/(dzeta^2);
LHS((i-1)*N+1, i*N+1) = L1(1,i)/(dx^2) + Lb2_LFL(i+1)*L3(1,i)/(4*dx*dzeta) + L5(1,i)/(2*dx);
end
end
for j = 1:N
LHS(j,j) = 1;
RHS(j,1) = 0;
LHS((Ms-1)*N+j,(Ms-1)*N+j) = 1;
RHS((Ms-1)*N+j,1) = 0;
end
v = LHS\RHS;
u_s = zeros(N, Ms);
for k = 1:N*Ms
if fix(k/N) == k/N
i = k/N;
else
i = fix(k/N) + 1;
end
j = k-(i-1)*N;
u_s(j,i) = v(k);
end
u_s(:, 1) = 0;
u_s(:, Ms) = 0;
end