-
Notifications
You must be signed in to change notification settings - Fork 3
/
coefficients_weno.m
109 lines (68 loc) · 1.55 KB
/
coefficients_weno.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
function [C]=coefficients_weno (r)
a=zeros(r,2);
C=zeros(r,2);
np=zeros(r,2);
nn=zeros(r,2);
%Derivada positiva (1)
x=0.5;
for k=0:r-1
for s=0:r
prod=1;
for l=0:r
if l~=s
prod=prod*((x)-(k-l+0.5));
end
end
a(k+1,1)=a(k+1,1)+prod;
end
if a(k+1,1)>0
np(1)=np(1)+1;
end
if a(k+1,1)<0
nn(1)=nn(1)+1;
end
end
for k=0:r-1
if a(k+1,1)==0
C(k+1,1)=1;
end
if a(k+1,1)>0
C(k+1,1)=1./(np(1)*a(k+1,1));
end
if a(k+1,1)<0
C(k+1,1)=1./(nn(1)*a(k+1,1));
end
end
%Derivada negativa (2)
x=-0.5;
for k=0:r-1
for s=0:r
prod=1;
for l=0:r
if l~=s
prod=prod*((x)-(k-l+0.5));
end
end
a(k+1,2)=a(k+1,2)+prod;
end
if a(k+1,2)>0
np(2)=np(2)+1;
end
if a(k+1,2)<0
nn(2)=nn(2)+1;
end
end
for k=0:r-1
if a(k+1,2)==0
C(k+1,2)=1;
end
if a(k+1,2)>0
C(k+1,2)=1./(np(2)*a(k+1,2));
end
if a(k+1,2)<0
C(k+1,2)=1./(nn(2)*a(k+1,2));
end
end
%**************************************************************************
C=abs(C);
end