2% MANUFACTUREDSOLUTION analytic density functions
3%
4% [sig,Nx0,h0,t0,mu] = manufacturedSolution(M,tFinal,s,beta)
5% returns in sig a cell array of M density function handles
6%
7% Input values:
8% M: number of spring scatterers (sources)
9% tFinal: final time
10% tol: error tolerance for the full problem (to compute starting dt)
11% gam: pad gam*Nyquist for window function (to compute starting dt)
12%
13% Outputs:
14% sig: function handle in terms of time that produces a vector of density
15% values evaluated at a specific time t:
16% [sig_1(t),...,sig_M(t)]
17% Nx0, h0: Number of grid points needed to ensure sig is resolved on the
18% time grid, as well as the maximum step-size.
19% t0, mu : (to document)
20% dataParam: is a struct involving t0,mu and sig.
21%
22% Notes:
23% Take the densities to be Gaussians for testing.
25override = 0; % turn on '1' to avoid restrictions on chosen test function
27mu = linspace(40,50,M)';
28t0 = linspace(1,3,M)';
30% choose t0 such that sig(0) = 0.01*eps
31% t0 = sqrt(log(amp/(0.01*eps))./mu);
33%%% Working resolution using standard deviation
34N0 = 20; % number of grid points per standard deviation (sd)
35bumpWidth = sqrt(log(1/eps)./mu);
36h0_vec = bumpWidth./N0;
37h0 = min(h0_vec);
39%%% Resolution computed from Fourier analysis
40% omega0 = 2*sqrt(mu).*sqrt(log((1/tol)*sqrt(pi./mu)));
41% h0 = min(pi*(1 - gam)./omega0);
42% Nt0 = ceil(tFinal/h0); % or Nt0 = ceil(2*pi/h0);
44% ensure that the forcing function is bandlimited within the computational time domain
45if(any(2*t0>tFinal)&&(override == 0))
46 error('Forcing inappropriate for testing. For tFinal=%1.2f\nChoose other mu or tFinal',tFinal);
47end
49% Define the density functions to be Gaussians
50sig = @(t) exp(-mu.*((t - t0).^2));
52dataParam.sig = sig;
53dataParam.t0 = t0;
54dataParam.mu = mu;
56end