Problem with finding roots of an interpolated function
4 ビュー (過去 30 日間)
古いコメントを表示
Alexander Engman
2016 年 10 月 17 日
編集済み: Alexander Engman
2016 年 10 月 17 日
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
2 件のコメント
Jan
2016 年 10 月 17 日
Splines are not necessarily a nice interpolation, because the trajectory can leave the range of the measurement data.
John D'Errico
2016 年 10 月 17 日
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.
採用された回答
その他の回答 (1 件)
John D'Errico
2016 年 10 月 17 日
編集済み: John D'Errico
2016 年 10 月 17 日
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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/175268/image.jpeg)
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)])
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/175269/image.jpeg)
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')
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/175270/image.jpeg)
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
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/175271/image.jpeg)
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:
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.
0 件のコメント
参考
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!