Thread Subject:

How do I solve an linear 1st order ode inverse problem?

Say I have the function: f(x) = x^2, 0 <= 1.
Instead of observing the function directly I am observing the derivative of the function + some noise (n): p(x) = df/dx + n.
Assume I also know that f(0) = 0.
Integrating the definition of p: I get: f(x) = int_0^x [p(y) - n]dy.
The problem with this approach is that noise also gets integrated as the script at the end of this post will show.

How can I estimate f(x) in a better way?

Thank in advance for any answers!

% inverseProblem.m
n = 100;
x = linspace(0,1,n);
dx = x(2) - x(1);

% "secret" function:
f = x.^2;

% observed derivative of f + random noise:
p = 2*x + randn(size(x));

subplot(1,2,1);
plot(x,2*x,'r');
hold on;
plot(x,p,'b');
xlabel('x');
legend('df/dx','p(x)');

subplot(1,2,2);
fi=cumsum(p)*dx;
plot(x,f,'r');
hold on
plot(x,fi,'b');
xlabel('x');
legend('f(x)','fi(x)');

Noise can be dealt with in several ways.
The first thing to ask is, do you know what type of noise is it? (shot noise, 1/f, etc..)
If the noise carries some distinct features you can try filter it out using functions from the signal processing toolbox, i.e. filter in the frequency domain etc..

Another generic approach is to smooth the signal (see the smooth function fir more details)
so you can try to use something like:
p = smooth(2*x + randn(size(x)),'sgolay');

Yet another approach is to use information regarding the function you are expecting to find. So the plot thickens as function of your function... :)

Just some generic though:

If you know some a priori information about the secret function, such as it belong to some functional set, e.g. polynomial, then you can use fitting methods.

Otherwise use different filtering. The mentioned Savitzki-Golay assumes the function to be locally polynomial.

Bruno