Many students ask me how do I do this or that in MATLAB. So I thought why not have a small series of my next few blogs do that. In this blog, I show you how to conduct spline interpolation.
- The MATLAB program link is here.
- The HTML version of the MATLAB program is here.
- DO NOT COPY AND PASTE THE PROGRAM BELOW BECAUSE THE SINGLE QUOTES DO NOT TRANSLATE TO THE CORRECT SINGLE QUOTES IN MATLAB EDITOR. DOWNLOAD THE MATLAB PROGRAM INSTEAD
%% HOW DO I DO THAT IN MATLAB SERIES?
% In this series, I am answering questions that students have asked
% me about MATLAB. Most of the questions relate to a mathematical
% procedure.
%% TOPIC
% How do I do spline interpolation?
%% SUMMARY
% Language : Matlab 2008a;
% Authors : Autar Kaw;
% Mfile available at
% http://numericalmethods.eng.usf.edu/blog/interpolate_spline.m;
% Last Revised : June 20, 2009;
% Abstract: This program shows you how to do spline interpolation?
% .
clc
clear all
clf
%% INTRODUCTION
disp(‘ABSTRACT’)
disp(‘ This program shows you how to do spline interpolation?’)
disp(‘ ‘)
disp(‘AUTHOR’)
disp(‘ Autar K Kaw of http://autarkaw.wordpress.com’)
disp(‘ ‘)
disp(‘MFILE SOURCE’)
disp(‘ http://numericalmethods.eng.usf.edu/blog/interpolation_spline.m’)
disp(‘ ‘)
disp(‘LAST REVISED’)
disp(‘ June 20, 2009′)
disp(‘ ‘)
%% INPUTS
% y vs x data to interpolate
% x data
x=[-1 -0.75 -0.5 -0.25 0.25 0.50 0.75 1];
% ydata
y=[-0.5 -0.5 -0.5 -0.5 0.5 0.5 0.5 0.5];
% Where do you want to interpolate at
xin=[-0.8 -0.7 0.7 0.8];
%% DISPLAYING INPUTS
disp(‘INPUTS’)
disp(‘The x data’)
x
disp(‘The y data’)
y
disp(‘The x values where you want to find the interpolated values’)
xin
disp(‘ ‘)
%% THE CODE
% Fitting to spline – it is cubic splines
yin=spline(x,y,xin);
% This is only for plotting the spline interpolants
% Find the number of data points
n=length(x);
xplot=x(1):(x(n)-x(1))/10000:x(n);
yplot=spline(x,y,xplot);
%% DISPLAYING OUTPUTS
disp(‘ ‘)
disp(‘OUTPUTS’)
disp(‘x values at which function is to be interpolated’)
xin
disp(‘y values at the xin values’)
yin
xlabel(‘x’);
ylabel(‘y’);
title(‘y vs x ‘);
plot(x,y,’o',’MarkerSize’,10,’MarkerEdgeColor’,'b’,'MarkerFaceColor’,'b’)
hold on
plot(xin,yin,’o',’MarkerSize’,10,’MarkerEdgeColor’,'r’,'MarkerFaceColor’,'r’)
hold on
plot(xplot,yplot,’LineWidth’,2)
legend(‘Points given’,'Points found’,'Spline Curve’,'Location’,'East’)
hold off
disp(‘ ‘)
This post is brought to you by Holistic Numerical Methods: Numerical Methods for the STEM undergraduate at http://numericalmethods.eng.usf.edu, the textbook on Numerical Methods with Applications available from the lulu storefront, and the YouTube video lectures available at http://numericalmethods.eng.usf.edu/videos and http://www.youtube.com/numericalmethodsguy
Subscribe to the blog via a reader or email to stay updated with this blog. Let the information follow you.
1 Comment
July 3, 2009 at 3:12 AM
Hi, my name is Mel. First of all, thank you very much in advance. And, I hope you could please response if you have any ideas/opinions regarding this mail as I’m a new Matlab user. So, I need your guidance, simple examples since you’re an expert.
Please, I would like to ask your help on how to code using Matlab for two – dimensional heat equation version based on the following one – dimensional heat equation version.
clear all; clc;
%n=number of points, Nstep=time-step
n=100;
Nstep=20;
%L=length of domain, k=heat diffusivity
L=1;
k=0.01;
%define dx and x
dx=2*L/(n-1);
dt=(0.46*dx^2)/k;
x=-L:dx:L;
%initial
for I=2:n-1
if I52
u(I)=0
else
u(I)=1;
end
end
u(1)=0;u(n)=0;
% store profile for computing and display
v=u;
v_initial=u;
title_str=’Heat equation:u_t-\kappa u_{xx}=0′;
%initialization u and v
plot(x,v_initial, ‘-’,'linewidth’,2);
set(gca, ‘fontsize’,14,’fontweight’,'bold’);
title(title_str);
axis tight;
set(gca,’nextplot’,'replacechildren’);
%time step
for time=1:Nstep
%end points
uxx(1)=v(2)-2*v(1)+v(n)
uxx(n)=v(1)-2*v(2)+v(n-1)
%interior points
for I=2:n-1
uxx(I)=v(I+1)-2*v(I)+v(I-1)
end
%rule
u=v+k*uxx;
%update v
v=u;
%plot
plot(x,v_initial,’-',x,u,’linewidth’,2);
picture(time)=getframe;
end
%display movie
disp(‘Playing animation…’);
movie(picture,1);