Problem with finding roots of an interpolated function

Hi!

I am having problem with an assignment. I am supposed to simulate a pendulum with Eulers method. I think I have managed to solve this part of the problem. The second step is to interpolate between the points with spline to create a function from which I can find the roots and thus find the period time of the pendulum (which would be the first root*2).

This is my code:

clear all, close all, clc

h=0.05; %step length

t=[0:h:2]; %interval length

g=9.82; %constant

l=0.5; %constant

x(1)=(25/180)*pi; %starting angle

y(1)=0; %starting velocity

T2=2*pi*sqrt(l/g) %actual period time

for i=1:length(t)-1

x(i+1)=x(i)+y(i)*h; %vector with angle values

y(i+1)=y(i)+(-g/l)*sin(x(i))*h; %vectore with velocity values

end;

f = @(t) spline(y,x,t); %interpolates the points in y and x

root=fzero(f,0.7)%finds roots of interpolated polynomial

T=root*2%period time

Tvec=[1.4182 1.1667 1.3815 1.3784 1.4055]; %vector with T-values from different h-values

for i=1:length(Tvec)

E2(i)=(abs(Tvec(i)-T2))/Tvec(i)%absolute error of period time

end;

plot(t,x,'-b'); %plots interpolated polynomial

My problem is that depending on the size of h, I am getting results of the root which vary from 0.7 (which roughly is accurate) to 0.55 (which is way off). I should be getting more and more accurate values as my h-value gets smaller and smaller. Since the error in the root does not seem to follow any kind of pattern, neither does the error of my period time. So in short, how should I rewrite this code in order to get correct values of the root of my interpolated function?

Your help will be much appreciated!

Thank you

Alexander Engman

[Jan]

Splines are not necessarily a nice interpolation, because the trajectory can leave the range of the measurement data.

[John D'Errico]

And this is why pchip MAY often be a better choice of interpolant. It works just like the function spline, but ensures the curve is bounded by the data while remaining fairly smooth.

Answer 1

[Alexander Engman]

Hi John!

Thank you for your very detailed and comprehensive answer.

I am afraid that I am not allowed to use the code that you have suggested. Is there any other way to make an interpolation for my polynomial that will actual work using any of the Matlab built in functions? My assignment clearly states that this is what I should do, so I guess there must be.

If I change my code to:

f = @(x2) spline(t,x,x2); %interpolates the points in t and x

ezplot(f,[min(y), max(y)])

Then I get something which actually resembles my pendulum and also includes the root I am looking for. However, fzero of this function at the approximate value of the root is giving me a completely different answer (around half the value of the actual root). Any recommendations on this issue?

Edit. In fact, just taking my value of the root *4 instead of *2 is giving me exactly the results I was hoping for. Is this a coincidence or is that right?

Thank you!

Alexander

Answer 2

[John D'Errico]

I think your problem is you do not understand spline interpolation, and what spline does when it sees your curve.

So I executed your code, up to the point where the spline is created. Before you go ANY further, always plot your data. LOOK AT WHAT YOU HAVE! THINK ABOUT WHAT IS HAPPENING!

...
plot(y,x,'o-')
grid on

You are using it in the function spline as if it is a function, x(y). This is true, since you defined f as

f = @(t) spline(y,x,t);

Note that in that plot, x is NOT a single valued function of y. Y is not a single valued function of x either, so the point is a bit moot. But spline does not know that, nor does it really care! It does what you ask it to do, and no more.

So now plot what f looks like:

ezplot(f,[min(y), max(y)])

UGH. See that this is NOT the original curve that you had. Spline took the points as ordered pairs (y(i),x(i)), and then tried to interpolate, but essentially sorting y. So it grabs each point in order of increasing y.

I think we can see what happened when we add the points to the spline plot.

hold on
plot(y,x,'ro')

Clearly you are not getting what you think you have.

The solution may be to use a tool like cscvn, which can deal with implicit function curves such as this. Or you could use my interparc tool, found on the file exchange. In any case, the issue is how to interpolate such a curve that is NOT a single valued function.

My recommendation (were this not surely homework) would be to use a tool like Doug Schwarz's intersections. It could find the intersections of your curve with another.

In this case, I'll use it where we have one curve that is your (x,y) pairs, in sequence from the original plot.

The second curve will be a straight line at x==0. So I'll connect the two points in the (x,y) plane: (0,-5) and (0,4).

[x0,y0] = intersections(x,y,[0 0],[-5 4])
Warning: NARGCHK will be removed in a future release. Use NARGINCHK or NARGOUTCHK instead. 
> In intersections (line 91) 
x0 =
     0
     0
     0
y0 =
         -4.36118590091745
         -2.26056827232397
          3.16741594889152

Sadly, Doug has not maintained this function, so it tosses out a warning message, but still works fine. There are three locations identified where the curve passes through x==0.

close
plot(y,x,'b-o')
grid on
hold on
plot(y0,x0,'rs-')
xlabel Y
ylabel X

As you can see, it has located the points of interest. Again, I've plotted them as you seem to be using them, with x as a function of y.

You can find Doug's very nice code here on the FEX:

http://www.mathworks.com/matlabcentral/fileexchange/11837-fast-and-robust-curve-intersections

Of course, if this is truly homework as I am sure it is, you will probably not be allowed to use Doug's code. But a simple linear interpolation will suffice.