1function [interpNodes,interpNodesIdx] = get_InterpNodes(auxNode, dt, order, endTime)
2% GETINTERPNODES retrieves the interpolation nodes on the time-grid to
3% interpolate unknown function values at auxiliary nodes
4%
5% [interpNodes,interpNodesIdx] = get_InterpNodes(auxNode, dt, order, endTime)
6% returns the interpolation nodes on the time-grid and their
7% corresponding indices to interpolate function values at a given
8% auxiliary node (off the time-grid). The function ensures that
9% interpolation nodes do not go beyong a given endTime
10%
11% Inputs:
12% auxNode: given value that may not coincide with the uniform time-grid
13% dt: time-step
14% order: equals the number of interpolation nodes used
15% endTime: ensure interpolation nodes < endTime
16%
17% Outputs:
18% interpNodes: a vector of size order that carries the interpolation nodes
19% needed to determine function values at auxNode
20% interpNodesIdx: a vector of size order that carries the interpolation
21% nodes indices on the time-grid
24%%% This is the one-sided treatment (towards the left)
25% rightGridPoint_idx = ceil(glnode/dt);
26% rightGridPoint = rightGridPoint_idx*dt; % find the nearest node on the uniform grid
27% interpNodes = linspace(rightGridPoint - (order - 1)*dt,rightGridPoint,order)';
28% interpNodesIdx = linspace((rightGridPoint_idx - (order - 1)),rightGridPoint_idx,order)';
30if nargin==0, test_getInterpNodes; return; end
32%%% This is the centered/skewed treatment
33m = order; % set m to equal the order
35tau_idx = round(auxNode/dt); % index of nearest node on time-grid
36tau = tau_idx*dt; % nearest node on time-grid
37endTime_idx = round(endTime/dt); % index of endTime
38k = endTime_idx - tau_idx; % how far is nearest node from endTime
40if(mod(m,2) == 0) % if m is even
41 if(k>=((m-2)/2)) % centered interpolation treatment
43 if(auxNode>tau)
44 % rightward skew: add interpolation point to the right end
45 startingNode = tau - ((m-2)/2)*dt;
46 startingIdx = tau_idx - ((m-2)/2);
48 % if last node is outside domain: add interpolation point to
49 % the left end
50 if((tau+((m/2)*dt))>endTime)
51 startingNode = tau - (m/2)*dt;
52 startingIdx = tau_idx - (m/2);
53 end
55 else % leftward skew: add interpolation point to the left end
56 startingNode = tau - (m/2)*dt;
57 startingIdx = tau_idx - (m/2);
58 end
59 else
60 % skewed treatment: when centered interpolation is not possible
61 % because interpolation nodes become larger than endTime
62 startingNode = tau - (m-k-1)*dt;
63 startingIdx = tau_idx - (m-k-1);
64 end
65else % if m is odd
66 if(k>=((m-1)/2)) % centered interpolation
67 startingNode = (tau - ((m-1)/2)*dt);
68 startingIdx = tau_idx - ((m-1)/2);
69 else % skewed treatment
70 startingNode = tau - (m-k-1)*dt;
71 startingIdx = tau_idx - (m-k-1);
72 end
73end
75% interpolation Nodes: startingNode, startingNode + dt,...
76% ...,startingNode + (m-1)*dt
78interpNodes = (startingNode:dt:(startingNode + ((m-1)*dt)))';
79interpNodesIdx = (startingIdx:(startingIdx + (m-1)))';
81end
83function test_getInterpNodes
84clf;
86startTime = 0; endTime = 1;
87t1 = .25; t2 = 1;
88N = 16;
89order = 5;
90saveToFile = 0; % set to '1'
91saveFile = '/Users/nalhassanieh/Desktop';
93tn = linspace(startTime,endTime,N);
94dt = tn(2) - tn(1);
95glnodes = lgwt(N,t1,t2);
97% plot integration domain and GL nodes
98figure(1);
99hAxes = axes('NextPlot','add',...
100 'DataAspectRatio',[1 1 1],...
101 'XLim',[startTime endTime],...
102 'YLim',[0 eps],...
103 'Color','none');
104p1 = plot(glnodes,0,'bo','MarkerSize',10);
105p2 = plot(tn,0,'r|','MarkerSize',10);
106p3 = plot([t1;t2],0,'kx','MarkerSize',10);
107l_temp = [p1(1);p2(1);p3(1)];
108% legend(l_temp,'GL nodes','time-grid','endpoints','interpreter','latex','location','north');
109% title('GL nodes for [t_1,t_2] on time-grid');
110set(gca,'XTick',[t1 t2], 'YTick', []);
111xticklabels({'$t_1$','$t_2$'})
112set(gca,'TickLabelInterpreter','latex');
114if(saveToFile == 1)
115 saveas(gcf,sprintf('%s/GLnodesOnTimeGrid',saveFile),'epsc');
116end
118pos = {'right','center','left'};
120cnt = 1;
121for glnodeNum = [2,8,15]
122 auxNode = glnodes(glnodeNum);
123 [interpNodes,~] = get_InterpNodes(auxNode, dt, order, endTime);
125 % plot one GL node and interpolation nodes
126 figure;
127 hAxes = axes('NextPlot','add',...
128 'DataAspectRatio',[1 1 1],...
129 'XLim',[startTime endTime],...
130 'YLim',[0 eps],...
131 'Color','none');
132 p1 = plot(glnodes(glnodeNum),0,'bo','MarkerSize',10);
133 p2 = plot(interpNodes,0,'*','color',[0.4660 0.6740 0.1880],'MarkerSize',10);
134 p3 = plot(tn,0,'r|','MarkerSize',10);
135 p4 = plot([t1;t2],0,'kx','MarkerSize',10);
136 l_temp = [p1(1);p2(1);p3(1);p4(1)];
137 % legend(l_temp,'GL node','interpolation nodes','time-grid','endpoints','interpreter','latex','location','north');
138 % title(sprintf('%s GL node interpolation',pos{cnt}));
139 set(gca,'XTick',[t1 t2], 'YTick', []);
140 xticklabels({'$t_1$','$t_2$'})
141 set(gca,'TickLabelInterpreter','latex');
143 if(saveToFile == 1)
144 saveas(gcf,sprintf('%s/exampleGLInterp_%s',saveFile,pos{cnt}),'epsc');
145 end
147 pause(0.1);
148 cnt = cnt + 1;
149end
151end