-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathDM_A1e3.m
68 lines (51 loc) · 1.78 KB
/
DM_A1e3.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
clc; clear all; close all;
k = 6; % Parameter in u(x)
x_max = 2*pi; x_min = 0;
x_an_step = 0.01;
x_an = x_min:x_an_step:x_max;
% Function handle for u(x) and its derivative
anSoln_handle = @(X, K) exp(K*sin(X));
du_exact_handle = @(X, K) K*cos(X) .* anSoln_handle(X, K);
% Threshold for minimum error
threshold = 1e-5;
N = 2; % Start with a small N
min_error = inf;
% Arrays to store N and min_error for plotting evolution
N_evolution = [];
min_error_evolution = [];
while min_error >= threshold
N = N + 2; % Increment N by 2 since it should be an even number
% odd method
j = 0:1:N;
x = (x_max/(N+1))*j;
[I, J] = ndgrid(j, j); % grid of indices
D = ((-1).^(I + J) ./ 2) .* (sin((I - J) .* pi ./ (N + 1))).^(-1);
D(1:(size(D,1)+1):end) = 0; % Set diagonal elements to zero
% Compute derivative using Fourier differentiation matrix
u = anSoln_handle(x, k);
du_approx = D * u';
du_real = du_exact_handle(x,k);
% Compute relative pointwise error
error = abs((du_real' - du_approx) ./ du_real');
% Measure minimum relative pointwise error
min_error = min(error);
% Store N and min_error for plotting evolution
N_evolution(end+1) = N;
min_error_evolution(end+1) = min_error;
% Plot
subplot(2,1,1)
plot(x_an, du_exact_handle(x_an,k), 'b', x, du_approx, 'r-');
xlabel('x');
ylabel('Derivative');
title(['Exact vs Approximate Derivative for N = ', num2str(N), ', k = ', num2str(k)]);
legend('Exact', 'Approximate');
drawnow;
subplot(2,1,2)
semilogy(N_evolution, min_error_evolution, 'r-');
hold on;
xlabel('N');
ylabel('Minimum Relative Error');
title(['Minimum Relative Error vs N for k = ', num2str(k)]);
drawnow;
end
fprintf('Minimum N value needed for error < 10^-5: %d\n', N);