1
$\begingroup$

I'm working on a problem in scientific computing namely fitting data to this equation

$c(z) = 4800 + p_1 + p_2 \cdot z/1000 + p_3 \cdot e^{ -p4 \cdot z/1000}$

The data is in a background question that I asked previously to prepared for the problem

acoustics under water

Now to split my work into smaller units, I should first try and fit the data. In the answer to the question is says that I can use lsqnonlin but it doesn't say how. I look in the manual for lsqnonlin and it says that the syntax is

x = lsqnonlin(fun,x0)

but I don't fully understand how to use it. Ideally I would read my data from a file, can you tell me how to do it? The data is

z c(z)
0 5050
500 4980
1000 4930
1500 4890
2000 4870
2500 4865
3000 4860
3500 4860
4000 4865
5000 4875
6000 4885
7000 4905
8000 4920
9000 4935
10000 4950
11000 4970
12000 4990

and the problem specification is problem U8 in this text.

Update

Running the code from the answer gives the following result for me after 400 iterations:

ans = 3.6277 14.8277 22.2631 -0.1114

So this means p1=3.6277, p2=14.8277, p3=22.2631, p4=-0.1114

Is that correct?

Update 2

>> z=0:10:12000;
>> fittedfun=4800 + 3.6277+14.8277/1000*z+22.2631*exp(0.1114/1000*z);
>> plot(z, fittedfun);

enter image description here

It looks rather linear, so I'm not totally certain that I'm seeing what I expect.

Update 3

Comparing the fitted function with the actual data the approximation doesn't look very good so it leads me to think that there is something wrong with my method or the equation.

enter image description here

>> plot(z, fittedfun);
>> z=0:500:12000;
>> fittedfun=4800 + 3.6277+14.8277/1000*z+22.2631*exp(0.1114/1000*z);
>> plot(z, fittedfun);
>> data=[5050 4980 4930 4890 4870 4865 4860 4860 4865 4875 4885 4905 4920 4935 4950 4970 4990];
>> z=[0:500:4000 5000:1000:12000];
>> zu=z

zu =

  Columns 1 through 6

           0         500        1000        1500        2000        2500

  Columns 7 through 12

        3000        3500        4000        5000        6000        7000

  Columns 13 through 17

        8000        9000       10000       11000       12000

>> plot(data,zu)
>> plot(zu,data)
>> hold on;
>> plot(zu,data)
>> fittedfun=4800 + 3.6277+14.8277/1000*z+22.2631*exp(0.1114/1000*z);
>> z=[0:500:4000 5000:1000:12000];
>> plot(z, fittedfun);
>> 
$\endgroup$
6
  • 1
    $\begingroup$ Yes, that is the correct answer (my MATLAB has the same output and the modeling is correct). Note, that I omited the (-) sign in the template function, so $p_4$ is to be read as $0.1114$. Whether the parameters are accurately fitted is up to you (I suggest plotting the function with those parameters and comparing the resulting curve to your expectations) $\endgroup$ – AlexR Jul 23 '13 at 12:51
  • 1
    $\begingroup$ I've found this article from the MathWorks on nonlinear fitting of exponential data helpful in the past. They use nlinfit rather than lsqnonlin, but I'd guess that the two are similar. $\endgroup$ – horchler Jul 23 '13 at 17:49
  • 1
    $\begingroup$ Your implementation still has the wrong sign in front of $p_4$. The function is supposed to read $$\ldots + 22.2631\cdot \exp(-0.1114 z)$$ $\endgroup$ – AlexR Aug 5 '13 at 20:20
  • 1
    $\begingroup$ @909 Niklas: is the Question still relevant to you? Please mark it as solved. $\endgroup$ – AlexR Aug 6 '13 at 13:24
  • 1
    $\begingroup$ @909Niklas I plotted the two on top of each other and it seemed okay to me. Also, MATLAB finished with "minimum found" and not with "maximal iterations reached", so I am confident, that this suffices. $\endgroup$ – AlexR Sep 2 '13 at 12:17
3
$\begingroup$

First, collapse the Function to only necessary parameters (all constants are "eaten up"): $$ c(z) = p_1 + p_2 z + p_3 \exp(p_4 z) = p_1 + p_2 z + p_3 \tilde{z}^{p_4}$$ Now iterate over $p_4$ in plausible values and do a 2D-polyfit with $x=z, y=\exp(z)^{p_4}$, see http://www.mathworks.com/matlabcentral/newsreader/view_thread/17164 (google matlab 2d polyfit). Finding $p_4$ can be tricky, though.

For the nonlinear solver, use an error measure (i.e. fun = @(pVec)norm(data-c(p1, p2, p3, p4)) where pVec(i)=p_i).

EDIT:
The following MATLAB snippet did the work for me, though used maxFunEvals (400) so no idea how good the approximation is:
z=[0:500:4000 5000:1000:12000];
data=[5050 4980 4930 4890 4870 4865 4860 4860 4865 4875 4885 4905 4920 4935 4950 4970 4990];
fun=@(p)4800 + p(1)+p(2)/1000*z+p(3)*exp(p(4)/1000*z)-data;
x0=[0 1 1 -1]; % Guess
lsqnonlin(fun,x0)
Another option would be to exploit $z=0$ to reduce the dimension if we force $c(0) = 5050$, but only if that is justified.

EDIT 2

The following script gives excellent results: z=[0:500:4000 5000:1000:12000]; data=[5050 4980 4930 4890 4870 4865 4860 4860 4865 4875 4885 4905 4920 4935 4950 4970 4990]; fun=@(p)(4800 + p(1))*ones(size(z)) +p(2)/1000*z+p(3)*exp(p(4)/1000*z)-data; x0=[0 0 0 -1]; % Guess opt = optimoptions('lsqnonlin', 'MaxFunEvals',1000); p=lsqnonlin(fun,x0,[],[],opt) fitf=@(t)(4800 + p(1))*ones(size(t)) + p(2)/1000*t+ p(3)*exp(p(4)/1000*t); tt=linspace(0,12000,1000); plot(z,data,'r-',tt,fitf(tt),'b-');

The only change was another guess in $x0$.

$\endgroup$
1
  • $\begingroup$ $e^{dz} \neq e^de^z = e^{d+z}.$ $\endgroup$ – Kirill Jul 23 '13 at 11:13

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.