diff --git a/lecture_17/octave-workspace b/lecture_17/octave-workspace index 52eed6b..c9ec2aa 100644 Binary files a/lecture_17/octave-workspace and b/lecture_17/octave-workspace differ diff --git a/lecture_18/Newtint.m b/lecture_18/Newtint.m index 93e77f6..e4c6c83 100644 --- a/lecture_18/Newtint.m +++ b/lecture_18/Newtint.m @@ -1,33 +1,34 @@ -function yint = Newtint_bak(x,y,xx) -% Newtint: Newton interpolating polynomial -% yint = Newtint(x,y,xx): Uses an (n - 1)-order Newton -% interpolating polynomial based on n data points (x, y) -% to determine a value of the dependent variable (yint) -% at a given value of the independent variable, xx. -% input: -% x = independent variable -% y = dependent variable -% xx = value of independent variable at which -% interpolation is calculated -% output: -% yint = interpolated value of dependent variable - -% compute the finite divided differences in the form of a -% difference table -n = length(x); -if length(y)~=n, error('x and y must be same length'); end -b = zeros(n,n); -% assign dependent variables to the first column of b. -b(:,1) = y(:); % the (:) ensures that y is a column vector. -for j = 2:n - for i = 1:n-j+1 - b(i,j) = (b(i+1,j-1)-b(i,j-1))/(x(i+j-1)-x(i)); - end -end -% use the finite divided differences to interpolate -xt = 1; -yint = b(1,1); -for j = 1:n-1 - xt = xt*(xx-x(j)); - yint = yint+b(1,j+1)*xt; -end +function yint = Newtint(x,y,xx) +% Newtint: Newton interpolating polynomial +% yint = Newtint(x,y,xx): Uses an (n - 1)-order Newton +% interpolating polynomial based on n data points (x, y) +% to determine a value of the dependent variable (yint) +% at a given value of the independent variable, xx. +% input: +% x = independent variable +% y = dependent variable +% xx = value of independent variable at which +% interpolation is calculated +% output: +% yint = interpolated value of dependent variable + +% compute the finite divided differences in the form of a +% difference table +n = length(x); +if length(y)~=n, error('x and y must be same length'); end +b = zeros(n,n); +% assign dependent variables to the first column of b. +b(:,1) = y(:); % the (:) ensures that y is a column vector. +for j = 2:n + for i = 1:n-j+1 + b(i,j) = (b(i+1,j-1)-b(i,j-1))/(x(i+j-1)-x(i)); + end +end +%b +% use the finite divided differences to interpolate +xt = 1; +yint = b(1,1); +for j = 1:n-1 + xt = xt*(xx-x(j)); + yint = yint+b(1,j+1)*xt; +end diff --git a/lecture_18/lecture_18.ipynb b/lecture_18/lecture_18.ipynb index 10a9d0a..bf7731c 100644 --- a/lecture_18/lecture_18.ipynb +++ b/lecture_18/lecture_18.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "code", - "execution_count": 1, + "execution_count": 146, "metadata": { "collapsed": true }, @@ -13,7 +13,7 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 107, "metadata": { "collapsed": true }, @@ -22,6 +22,40 @@ "%plot --format svg" ] }, + { + "cell_type": "code", + "execution_count": 115, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ans =\n", + "\n", + " 0.081014\n", + " 0.317493\n", + " 0.690279\n", + " 1.169170\n", + " 1.715370\n", + " 2.284630\n", + " 2.830830\n", + " 3.309721\n", + " 3.682507\n", + " 3.918986\n", + "\n" + ] + } + ], + "source": [ + "N=10;\n", + "A_beam=diag(ones(N,1))*2+diag(ones(N-1,1)*-1,-1)+diag(ones(N-1,1)*-1,1);\n", + "[v,e]=eig(A_beam);\n", + "diag(e)" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -108,7 +142,7 @@ }, { "cell_type": "code", - "execution_count": 16, + "execution_count": 116, "metadata": { "collapsed": false }, @@ -293,7 +327,7 @@ }, { "cell_type": "code", - "execution_count": 18, + "execution_count": 127, "metadata": { "collapsed": false }, @@ -316,7 +350,7 @@ }, { "cell_type": "code", - "execution_count": 19, + "execution_count": 128, "metadata": { "collapsed": false }, @@ -490,6 +524,167 @@ "plot(data(:,1),data(:,2),'o',data(:,1),yhat)" ] }, + { + "cell_type": "code", + "execution_count": 130, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t-15\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-10\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t15\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t6\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t8\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t10\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\t\n", + "\t\tresiduals of function\n", + "\t\n", + "\n", + "\n", + "\n", + "\tgnuplot_plot_1a\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "[sse,yhat]=sse_nonlin_exp(a,data(:,1),data(:,2));\n", + "plot(data(:,1),data(:,2)-yhat)\n", + "title('residuals of function')" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -506,12 +701,12 @@ "\n", "$\\sigma(t)=\\frac{1}{1+e^{-t}}$\n", "\n", - "We can use this function to describe the likelihood of failure (1) or success (0). When z=0, the probability of failure is 50%. " + "We can use this function to describe the likelihood of failure (1) or success (0). When t=0, the probability of failure is 50%. " ] }, { "cell_type": "code", - "execution_count": 30, + "execution_count": 131, "metadata": { "collapsed": false }, @@ -691,7 +886,7 @@ }, { "cell_type": "code", - "execution_count": 25, + "execution_count": 132, "metadata": { "collapsed": false }, @@ -739,7 +934,7 @@ }, { "cell_type": "code", - "execution_count": 27, + "execution_count": 134, "metadata": { "collapsed": false }, @@ -786,115 +981,125 @@ "\n", "\n", "\t\n", - "\t\t\n", + "\t\t\n", "\t\n", "\n", "\n", "\n", "\n", - "\t\t\n", + "\t\t\n", "\t\t0\n", "\t\n", "\n", "\n", - "\t\t\n", + "\t\t\n", "\t\t0.2\n", "\t\n", "\n", "\n", - "\t\t\n", + "\t\t\n", "\t\t0.4\n", "\t\n", "\n", "\n", - "\t\t\n", + "\t\t\n", "\t\t0.6\n", "\t\n", "\n", "\n", - "\t\t\n", + "\t\t\n", "\t\t0.8\n", "\t\n", "\n", "\n", - "\t\t\n", + "\t\t\n", "\t\t1\n", "\t\n", "\n", "\n", - "\t\t\n", + "\t\t\n", "\t\t50\n", "\t\n", "\n", "\n", - "\t\t\n", + "\t\t\n", "\t\t55\n", "\t\n", "\n", "\n", - "\t\t\n", + "\t\t\n", "\t\t60\n", "\t\n", "\n", "\n", - "\t\t\n", + "\t\t\n", "\t\t65\n", "\t\n", "\n", "\n", - "\t\t\n", + "\t\t\n", "\t\t70\n", "\t\n", "\n", "\n", - "\t\t\n", + "\t\t\n", "\t\t75\n", "\t\n", "\n", "\n", - "\t\t\n", + "\t\t\n", "\t\t80\n", "\t\n", "\n", "\n", - "\t\t\n", + "\t\t\n", "\t\t85\n", "\t\n", "\n", "\n", "\n", "\n", - "\t\n", + "\t\n", + "\n", + "\t\n", + "\t\tfailure (1)/ pass (0)\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\tTemp (F)\n", + "\t\n", + "\n", "\n", "\n", "\tgnuplot_plot_1a\n", - "\n", + "\n", "\n", "\n", "\t \n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", "\n", "\t\n", "\n", @@ -917,19 +1122,21 @@ } ], "source": [ - "plot(oring(:,2),oring(:,3),'o')" + "plot(oring(:,2),oring(:,3),'o')\n", + "xlabel('Temp (F)')\n", + "ylabel('failure (1)/ pass (0)')" ] }, { "cell_type": "code", - "execution_count": 54, + "execution_count": 135, "metadata": { "collapsed": true }, "outputs": [], "source": [ "function J=sse_logistic(a,x,y)\n", - " % Create function to calculate SSE of logistic function\n", + " % Create function to calculate cost of logistic function\n", " % t = a0+a1*x\n", " % sigma(t) = 1./(1+e^(-t))\n", " sigma=@(t) 1./(1+exp(-t));\n", @@ -942,7 +1149,7 @@ }, { "cell_type": "code", - "execution_count": 61, + "execution_count": 142, "metadata": { "collapsed": false }, @@ -1006,36 +1213,41 @@ "\n", "\n", "\n", - "\t\t\n", + "\t\t\n", "\t\t0\n", "\t\n", "\n", "\n", - "\t\t\n", + "\t\t\n", "\t\t0.2\n", "\t\n", "\n", "\n", - "\t\t\n", + "\t\t\n", "\t\t0.4\n", "\t\n", "\n", "\n", - "\t\t\n", + "\t\t\n", "\t\t0.6\n", "\t\n", "\n", "\n", - "\t\t\n", + "\t\t\n", "\t\t0.8\n", "\t\n", "\n", "\n", - "\t\t\n", + "\t\t\n", "\t\t1\n", "\t\n", "\n", "\n", + "\t\t\n", + "\t\t1.2\n", + "\t\n", + "\n", + "\n", "\t\t\n", "\t\t50\n", "\t\n", @@ -1086,36 +1298,40 @@ "\n", "\n", "\t \n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", "\n", "\t\n", "\tgnuplot_plot_2a\n", "\n", - "\t\n", + "\t\n", "\t\n", - "\n", + "\tgnuplot_plot_3a\n", + "\n", + "\t\n", + "\t\n", + "\n", "\n", "\n", "\n", @@ -1139,12 +1355,13 @@ "a=fminsearch(@(a) sse_logistic(a,oring(:,2),oring(:,3)),[0,-3])\n", "\n", "T=linspace(50,85);\n", - "plot(oring(:,2),oring(:,3),'o',T,sigma(a(1)+a(2)*T))" + "plot(oring(:,2),oring(:,3),'o',T,sigma(a(1)+a(2)*T),T,a(1)+a(2)*T)\n", + "axis([50,85,-0.1,1.2])" ] }, { "cell_type": "code", - "execution_count": 75, + "execution_count": 139, "metadata": { "collapsed": false }, @@ -1154,13 +1371,15 @@ "output_type": "stream", "text": [ "probability of failure when 70 degrees is 23.00% \n", - "probability of failure when 60 degrees is 75.25%\n" + "probability of failure when 60 degrees is 75.25%\n", + "probability of failure when 36 degrees is 99.87%\n" ] } ], "source": [ "fprintf('probability of failure when 70 degrees is %1.2f%% ',100*sigma(a(1)+a(2)*70))\n", - "fprintf('probability of failure when 60 degrees is %1.2f%%',100*sigma(a(1)+a(2)*60))\n" + "fprintf('probability of failure when 60 degrees is %1.2f%%',100*sigma(a(1)+a(2)*60))\n", + "fprintf('probability of failure when 36 degrees is %1.2f%%',100*sigma(a(1)+a(2)*36))\n" ] }, { @@ -1196,7 +1415,7 @@ }, { "cell_type": "code", - "execution_count": 84, + "execution_count": 149, "metadata": { "collapsed": false }, @@ -1207,6 +1426,7 @@ "text": [ "ln(2)~0.358352\n", "ln(2)~0.462098\n", + "ln(2)~0.549306\n", "ln(2)=0.693147\n" ] } @@ -1215,14 +1435,16 @@ "ln2_16=log(1)+(log(6)-log(1))/(6-1)*(2-1);\n", "fprintf('ln(2)~%f\\n',ln2_16)\n", "ln2_14=log(1)+(log(4)-log(1))/(4-1)*(2-1);\n", + "ln2_13=log(1)+(log(3)-log(1))/(3-1)*(2-1);\n", "fprintf('ln(2)~%f\\n',ln2_14)\n", + "fprintf('ln(2)~%f\\n',ln2_13)\n", "ln2=log(2);\n", "fprintf('ln(2)=%f\\n',ln2)" ] }, { "cell_type": "code", - "execution_count": 87, + "execution_count": 147, "metadata": { "collapsed": false }, @@ -1269,98 +1491,108 @@ "\n", "\n", "\t\n", - "\t\t\n", + "\t\t\n", "\t\n", "\n", "\n", "\n", "\n", - "\t\t\n", + "\t\t\n", "\t\t0\n", "\t\n", "\n", "\n", - "\t\t\n", + "\t\t\n", "\t\t0.5\n", "\t\n", "\n", "\n", - "\t\t\n", + "\t\t\n", "\t\t1\n", "\t\n", "\n", "\n", - "\t\t\n", + "\t\t\n", "\t\t1.5\n", "\t\n", "\n", "\n", - "\t\t\n", + "\t\t\n", "\t\t2\n", "\t\n", "\n", "\n", - "\t\t\n", + "\t\t\n", "\t\t1\n", "\t\n", "\n", "\n", - "\t\t\n", + "\t\t\n", "\t\t2\n", "\t\n", "\n", "\n", - "\t\t\n", + "\t\t\n", "\t\t3\n", "\t\n", "\n", "\n", - "\t\t\n", + "\t\t\n", "\t\t4\n", "\t\n", "\n", "\n", - "\t\t\n", + "\t\t\n", "\t\t5\n", "\t\n", "\n", "\n", - "\t\t\n", + "\t\t\n", "\t\t6\n", "\t\n", "\n", "\n", "\n", "\n", - "\t\n", + "\t\n", + "\n", + "\t\n", + "\t\tln(x)\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\tx\n", + "\t\n", + "\n", "\n", "\n", "\tgnuplot_plot_1a\n", - "\n", + "\n", "\n", "\n", - "\t\n", + "\t\n", "\t\n", "\tgnuplot_plot_2a\n", "\n", "\t \n", - "\t\n", + "\t\n", "\n", "\t\n", "\tgnuplot_plot_3a\n", "\n", - "\t\t \n", - "\t\n", - "\t\n", - "\t\n", + "\t\t \n", + "\t\n", + "\t\n", + "\t\n", "\n", "\t\n", "\tgnuplot_plot_4a\n", "\n", - "\t\t \n", - "\t\n", - "\t\n", - "\t\n", + "\t\t \n", + "\t\n", + "\t\n", + "\t\n", "\n", "\t\n", "\n", @@ -1386,7 +1618,9 @@ "x=linspace(1,6);\n", "plot(x,log(x),2,log(2),'*',...\n", "[1,2,6],[log(1),ln2_16,log(6)],'o-',...\n", - "[1,2,4],[log(1),ln2_14,log(4)],'s-')" + "[1,2,4],[log(1),ln2_14,log(4)],'s-')\n", + "ylabel('ln(x)')\n", + "xlabel('x')" ] }, { @@ -1417,7 +1651,41 @@ }, { "cell_type": "code", - "execution_count": 89, + "execution_count": 154, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Z =\n", + "\n", + " 1 1 1\n", + " 1 4 16\n", + " 1 600 360000\n", + "\n", + "ans = 5.1766e+05\n", + "ans =\n", + "\n", + " -4.6513e-01\n", + " 4.6589e-01\n", + " -7.5741e-04\n", + "\n" + ] + } + ], + "source": [ + "x=[1,4,600]';\n", + "Z=[x.^0,x.^1,x.^2]\n", + "cond(Z)\n", + "Z\\log(x)" + ] + }, + { + "cell_type": "code", + "execution_count": 155, "metadata": { "collapsed": false }, @@ -1441,18 +1709,25 @@ "f3=log(x3);\n", "\n", "b1=f1\n", - "b2=(f2-f1)/(x2-x1)\n", + "b2=(f2-b1)/(x2-x1)\n", "b3=(f3-f2)/(x3-x2)-b2;\n", "b3=b3/(x3-x1)" ] }, { "cell_type": "code", - "execution_count": 91, + "execution_count": 157, "metadata": { "collapsed": false }, "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ans = 0.56584\r\n" + ] + }, { "data": { "image/svg+xml": [ @@ -1615,7 +1890,8 @@ "f=@(x) b1+b2*(x-x1)+b3*(x-x1).*(x-x2);\n", "plot(x,log(x),2,log(2),'*',...\n", "[1,4,6],[log(1),log(4),log(6)],'ro',...\n", - "x,f(x),'r-',2,f(2),'s')" + "x,f(x),'r-',2,f(2),'s')\n", + "f(2)" ] }, { @@ -1649,7 +1925,41 @@ }, { "cell_type": "code", - "execution_count": 105, + "execution_count": 160, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "b =\n", + "\n", + " 0.00000 0.54931 -0.08721 0.01178\n", + " 1.09861 0.28768 -0.02832 0.00000\n", + " 1.38629 0.20273 0.00000 0.00000\n", + " 1.79176 0.00000 0.00000 0.00000\n", + "\n", + "ans = 0.66007\n", + "ans =\n", + "\n", + " 0.00000\n", + " 1.09861\n", + " 1.38629\n", + " 1.79176\n", + "\n" + ] + } + ], + "source": [ + "Newtint([1,3,4,6],log([1,3,4,6]),2)\n", + "log([1,3,4,6]')" + ] + }, + { + "cell_type": "code", + "execution_count": 162, "metadata": { "collapsed": false }, @@ -1659,7 +1969,7 @@ "output_type": "stream", "text": [ "ln(2)=0.693147\n", - "ln(2)~0.722462\n" + "ln(2)~0.366349\n" ] }, { @@ -1704,113 +2014,113 @@ "\n", "\n", "\t\n", - "\t\t\n", + "\t\t\n", "\t\n", "\n", "\n", "\n", "\n", - "\t\t\n", - "\t\t-1\n", + "\t\t\n", + "\t\t-2\n", "\t\n", "\n", "\n", - "\t\t\n", - "\t\t-0.5\n", + "\t\t\n", + "\t\t0\n", "\t\n", "\n", "\n", - "\t\t\n", - "\t\t0\n", + "\t\t\n", + "\t\t2\n", "\t\n", "\n", "\n", - "\t\t\n", - "\t\t0.5\n", + "\t\t\n", + "\t\t4\n", "\t\n", "\n", "\n", - "\t\t\n", - "\t\t1\n", + "\t\t\n", + "\t\t6\n", "\t\n", "\n", "\n", - "\t\t\n", - "\t\t1.5\n", + "\t\t\n", + "\t\t8\n", "\t\n", "\n", "\n", - "\t\t\n", - "\t\t2\n", + "\t\t\n", + "\t\t10\n", "\t\n", "\n", "\n", - "\t\t\n", + "\t\t\n", "\t\t0\n", "\t\n", "\n", "\n", - "\t\t\n", - "\t\t1\n", + "\t\t\n", + "\t\t10\n", "\t\n", "\n", "\n", - "\t\t\n", - "\t\t2\n", + "\t\t\n", + "\t\t20\n", "\t\n", "\n", "\n", - "\t\t\n", - "\t\t3\n", + "\t\t\n", + "\t\t30\n", "\t\n", "\n", "\n", - "\t\t\n", - "\t\t4\n", + "\t\t\n", + "\t\t40\n", "\t\n", "\n", "\n", - "\t\t\n", - "\t\t5\n", + "\t\t\n", + "\t\t50\n", "\t\n", "\n", "\n", "\t\t\n", - "\t\t6\n", + "\t\t60\n", "\t\n", "\n", "\n", "\n", "\n", - "\t\n", + "\t\n", "\n", "\n", "\tgnuplot_plot_1a\n", "\n", "\n", "\n", - "\t\n", + "\t\n", "\t\n", "\tgnuplot_plot_2a\n", "\n", "\t \n", - "\t\n", + "\t\n", "\n", "\t\n", "\tgnuplot_plot_3a\n", "\n", "\t \n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", - "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", "\n", "\t\n", "\tgnuplot_plot_4a\n", "\n", - "\t\n", + "\t\n", "\t\n", "\n", "\n", @@ -1833,10 +2143,11 @@ ], "source": [ "\n", - "yy=zeros(size(xx));\n", - "x=[0.5,1,3,4,5,6]; % define independent var's\n", + "\n", + "x=[0.2,3,10,20,50,60]; % define independent var's\n", "y=log(x); % define dependent var's\n", "xx=linspace(min(x),max(x));\n", + "yy=zeros(size(xx));\n", "for i=1:length(xx)\n", " yy(i)=Newtint(x,y,xx(i));\n", "end\n", diff --git a/lecture_18/octave-workspace b/lecture_18/octave-workspace new file mode 100644 index 0000000..c651696 Binary files /dev/null and b/lecture_18/octave-workspace differ diff --git a/lecture_19/.ipynb_checkpoints/lecture 19-checkpoint.ipynb b/lecture_19/.ipynb_checkpoints/lecture 19-checkpoint.ipynb new file mode 100644 index 0000000..2911efe --- /dev/null +++ b/lecture_19/.ipynb_checkpoints/lecture 19-checkpoint.ipynb @@ -0,0 +1,497 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "setdefaults" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "%plot --format svg" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Splines (Brief introduction before next section)\n", + "\n", + "Following interpolation discussion, instead of estimating 9 data points with an eighth-order polynomial, it makes more sense to fit sections of the curve to lower-order polynomials:\n", + "\n", + "0. zeroth-order (nearest neighbor)\n", + "\n", + "1. first-order (linear interpolation)\n", + "\n", + "2. third-order (cubic interpolation)\n", + "\n", + "Matlab and Octave have built-in functions for 1D and 2D interpolation:\n", + "\n", + "`interp1`\n", + "\n", + "`interp2`" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "'interp1' is a function from the file /usr/share/octave/4.0.0/m/general/interp1.m\n", + "\n", + " -- Function File: YI = interp1 (X, Y, XI)\n", + " -- Function File: YI = interp1 (Y, XI)\n", + " -- Function File: YI = interp1 (..., METHOD)\n", + " -- Function File: YI = interp1 (..., EXTRAP)\n", + " -- Function File: YI = interp1 (..., \"left\")\n", + " -- Function File: YI = interp1 (..., \"right\")\n", + " -- Function File: PP = interp1 (..., \"pp\")\n", + "\n", + " One-dimensional interpolation.\n", + "\n", + " Interpolate input data to determine the value of YI at the points\n", + " XI. If not specified, X is taken to be the indices of Y ('1:length\n", + " (Y)'). If Y is a matrix or an N-dimensional array, the\n", + " interpolation is performed on each column of Y.\n", + "\n", + " The interpolation METHOD is one of:\n", + "\n", + " \"nearest\"\n", + " Return the nearest neighbor.\n", + "\n", + " \"previous\"\n", + " Return the previous neighbor.\n", + "\n", + " \"next\"\n", + " Return the next neighbor.\n", + "\n", + " \"linear\" (default)\n", + " Linear interpolation from nearest neighbors.\n", + "\n", + " \"pchip\"\n", + " Piecewise cubic Hermite interpolating\n", + " polynomial--shape-preserving interpolation with smooth first\n", + " derivative.\n", + "\n", + " \"cubic\"\n", + " Cubic interpolation (same as \"pchip\").\n", + "\n", + " \"spline\"\n", + " Cubic spline interpolation--smooth first and second\n", + " derivatives throughout the curve.\n", + "\n", + " Adding '*' to the start of any method above forces 'interp1' to\n", + " assume that X is uniformly spaced, and only 'X(1)' and 'X(2)' are\n", + " referenced. This is usually faster, and is never slower. The\n", + " default method is \"linear\".\n", + "\n", + " If EXTRAP is the string \"extrap\", then extrapolate values beyond\n", + " the endpoints using the current METHOD. If EXTRAP is a number,\n", + " then replace values beyond the endpoints with that number. When\n", + " unspecified, EXTRAP defaults to 'NA'.\n", + "\n", + " If the string argument \"pp\" is specified, then XI should not be\n", + " supplied and 'interp1' returns a piecewise polynomial object. This\n", + " object can later be used with 'ppval' to evaluate the\n", + " interpolation. There is an equivalence, such that 'ppval (interp1\n", + " (X, Y, METHOD, \"pp\"), XI) == interp1 (X, Y, XI, METHOD, \"extrap\")'.\n", + "\n", + " Duplicate points in X specify a discontinuous interpolant. There\n", + " may be at most 2 consecutive points with the same value. If X is\n", + " increasing, the default discontinuous interpolant is\n", + " right-continuous. If X is decreasing, the default discontinuous\n", + " interpolant is left-continuous. The continuity condition of the\n", + " interpolant may be specified by using the options \"left\" or \"right\"\n", + " to select a left-continuous or right-continuous interpolant,\n", + " respectively. Discontinuous interpolation is only allowed for\n", + " \"nearest\" and \"linear\" methods; in all other cases, the X-values\n", + " must be unique.\n", + "\n", + " An example of the use of 'interp1' is\n", + "\n", + " xf = [0:0.05:10];\n", + " yf = sin (2*pi*xf/5);\n", + " xp = [0:10];\n", + " yp = sin (2*pi*xp/5);\n", + " lin = interp1 (xp, yp, xf);\n", + " near = interp1 (xp, yp, xf, \"nearest\");\n", + " pch = interp1 (xp, yp, xf, \"pchip\");\n", + " spl = interp1 (xp, yp, xf, \"spline\");\n", + " plot (xf,yf,\"r\", xf,near,\"g\", xf,lin,\"b\", xf,pch,\"c\", xf,spl,\"m\",\n", + " xp,yp,\"r*\");\n", + " legend (\"original\", \"nearest\", \"linear\", \"pchip\", \"spline\");\n", + "\n", + " See also: pchip, spline, interpft, interp2, interp3, interpn.\n", + "\n", + "Additional help for built-in functions and operators is\n", + "available in the online version of the manual. Use the command\n", + "'doc ' to search the manual index.\n", + "\n", + "Help and information about Octave is also available on the WWW\n", + "at http://www.octave.org and via the help@octave.org\n", + "mailing list.\n" + ] + } + ], + "source": [ + "help interp1" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "'interp2' is a function from the file /usr/share/octave/4.0.0/m/general/interp2.m\n", + "\n", + " -- Function File: ZI = interp2 (X, Y, Z, XI, YI)\n", + " -- Function File: ZI = interp2 (Z, XI, YI)\n", + " -- Function File: ZI = interp2 (Z, N)\n", + " -- Function File: ZI = interp2 (Z)\n", + " -- Function File: ZI = interp2 (..., METHOD)\n", + " -- Function File: ZI = interp2 (..., METHOD, EXTRAP)\n", + "\n", + " Two-dimensional interpolation.\n", + "\n", + " Interpolate reference data X, Y, Z to determine ZI at the\n", + " coordinates XI, YI. The reference data X, Y can be matrices, as\n", + " returned by 'meshgrid', in which case the sizes of X, Y, and Z must\n", + " be equal. If X, Y are vectors describing a grid then 'length (X)\n", + " == columns (Z)' and 'length (Y) == rows (Z)'. In either case the\n", + " input data must be strictly monotonic.\n", + "\n", + " If called without X, Y, and just a single reference data matrix Z,\n", + " the 2-D region 'X = 1:columns (Z), Y = 1:rows (Z)' is assumed.\n", + " This saves memory if the grid is regular and the distance between\n", + " points is not important.\n", + "\n", + " If called with a single reference data matrix Z and a refinement\n", + " value N, then perform interpolation over a grid where each original\n", + " interval has been recursively subdivided N times. This results in\n", + " '2^N-1' additional points for every interval in the original grid.\n", + " If N is omitted a value of 1 is used. As an example, the interval\n", + " [0,1] with 'N==2' results in a refined interval with points at [0,\n", + " 1/4, 1/2, 3/4, 1].\n", + "\n", + " The interpolation METHOD is one of:\n", + "\n", + " \"nearest\"\n", + " Return the nearest neighbor.\n", + "\n", + " \"linear\" (default)\n", + " Linear interpolation from nearest neighbors.\n", + "\n", + " \"pchip\"\n", + " Piecewise cubic Hermite interpolating\n", + " polynomial--shape-preserving interpolation with smooth first\n", + " derivative.\n", + "\n", + " \"cubic\"\n", + " Cubic interpolation (same as \"pchip\").\n", + "\n", + " \"spline\"\n", + " Cubic spline interpolation--smooth first and second\n", + " derivatives throughout the curve.\n", + "\n", + " EXTRAP is a scalar number. It replaces values beyond the endpoints\n", + " with EXTRAP. Note that if EXTRAPVAL is used, METHOD must be\n", + " specified as well. If EXTRAP is omitted and the METHOD is\n", + " \"spline\", then the extrapolated values of the \"spline\" are used.\n", + " Otherwise the default EXTRAP value for any other METHOD is \"NA\".\n", + "\n", + " See also: interp1, interp3, interpn, meshgrid.\n", + "\n", + "Additional help for built-in functions and operators is\n", + "available in the online version of the manual. Use the command\n", + "'doc ' to search the manual index.\n", + "\n", + "Help and information about Octave is also available on the WWW\n", + "at http://www.octave.org and via the help@octave.org\n", + "mailing list.\n" + ] + } + ], + "source": [ + "help interp2" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t-1.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-0.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-3\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t3\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\tdata\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\tdata\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\n", + "\t\n", + "\tlinear\n", + "\n", + "\t\n", + "\t\tlinear\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tcubic spline\n", + "\n", + "\t\n", + "\t\tcubic spline\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tpiecewise cubic\n", + "\n", + "\t\n", + "\t\tpiecewise cubic\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "x=linspace(-pi,pi,9);\n", + "xi=linspace(-pi,pi,100);\n", + "y=sin(x);\n", + "yi_lin=interp1(x,y,xi,'linear');\n", + "yi_spline=interp1(x,y,xi,'spline'); \n", + "yi_cubic=interp1(x,y,xi,'cubic');\n", + "plot(x,y,'o',xi,yi_lin,xi,yi_spline,xi,yi_cubic)\n", + "axis([-pi,pi,-1.5,1.5])\n", + "legend('data','linear','cubic spline','piecewise cubic','Location','NorthWest')\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Example: Accelerate then hold velocity\n", + "\n", + "Here the time is given as vector t in seconds and the " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "t=[0 2 40 56 68 80 84 96 104 110]';\n", + "v=[0 20 20 38 80 80 100 100 125 125]';\n" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Octave", + "language": "octave", + "name": "octave" + }, + "language_info": { + "file_extension": ".m", + "help_links": [ + { + "text": "MetaKernel Magics", + "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md" + } + ], + "mimetype": "text/x-octave", + "name": "octave", + "version": "0.19.14" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/lecture_19/Newtint.m b/lecture_19/Newtint.m new file mode 100644 index 0000000..e4c6c83 --- /dev/null +++ b/lecture_19/Newtint.m @@ -0,0 +1,34 @@ +function yint = Newtint(x,y,xx) +% Newtint: Newton interpolating polynomial +% yint = Newtint(x,y,xx): Uses an (n - 1)-order Newton +% interpolating polynomial based on n data points (x, y) +% to determine a value of the dependent variable (yint) +% at a given value of the independent variable, xx. +% input: +% x = independent variable +% y = dependent variable +% xx = value of independent variable at which +% interpolation is calculated +% output: +% yint = interpolated value of dependent variable + +% compute the finite divided differences in the form of a +% difference table +n = length(x); +if length(y)~=n, error('x and y must be same length'); end +b = zeros(n,n); +% assign dependent variables to the first column of b. +b(:,1) = y(:); % the (:) ensures that y is a column vector. +for j = 2:n + for i = 1:n-j+1 + b(i,j) = (b(i+1,j-1)-b(i,j-1))/(x(i+j-1)-x(i)); + end +end +%b +% use the finite divided differences to interpolate +xt = 1; +yint = b(1,1); +for j = 1:n-1 + xt = xt*(xx-x(j)); + yint = yint+b(1,j+1)*xt; +end diff --git a/lecture_19/lecture 19.ipynb b/lecture_19/lecture 19.ipynb new file mode 100644 index 0000000..e6f9f76 --- /dev/null +++ b/lecture_19/lecture 19.ipynb @@ -0,0 +1,1488 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "setdefaults" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "%plot --format svg" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Questions from last class\n", + "\n", + "![q1](q1.png)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "![q2](q2.png)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "![q3](q3.png)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "![q4](q4.png)\n", + "\n", + "#### Other:\n", + "\n", + "Twice the amount of points needed\n", + "\n", + "depends on what order polynomial it is and how far the data needs to be extrapolated \n", + "\n", + "As man you as possible \n", + "\n", + "Never extrapolate unless linear interpolation.\n", + "\n", + "You shouldn't. 2 if the linear is a good fit for the region, and you absolutely have to.\n", + "\n", + "Wait can you do that?\n", + "\n", + "Don't use extrapolation\n", + "\n", + "do not extrapolate\n", + "\n", + "As many data points as you have\n", + "\n", + "the more the better so that the best polynomial can be made through the data points\n", + "\n", + "Twice the amount of points needed" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Questions from you\n", + "\n", + "- when will the project assignment be finalized? Also do you pronounce it \"jiff\" or \"gif\"?\n", + "\n", + "- If blue is red and red is blue, then what is purple? \n", + "\n", + "- How do we open the .ipynb lecture files? Or will the lectures continue to be also saved in pdf (last few have not).\n", + "\n", + "- When will we be put on teams for the final project?\n", + "\n", + "- What is the grading rubric for the project?\n", + "\n", + "- How to sync repository with files from laptop like hw without using Github desktop \n", + "\n", + "- Are there any upcoming deadlines for the project?\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Splines (Brief introduction before next section)\n", + "\n", + "Following interpolation discussion, instead of estimating 9 data points with an eighth-order polynomial, it makes more sense to fit sections of the curve to lower-order polynomials:\n", + "\n", + "0. zeroth-order (nearest neighbor)\n", + "\n", + "1. first-order (linear interpolation)\n", + "\n", + "2. third-order (cubic interpolation)\n", + "\n", + "Matlab and Octave have built-in functions for 1D and 2D interpolation:\n", + "\n", + "`interp1`\n", + "\n", + "`interp2`" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "'interp1' is a function from the file /usr/share/octave/4.0.0/m/general/interp1.m\n", + "\n", + " -- Function File: YI = interp1 (X, Y, XI)\n", + " -- Function File: YI = interp1 (Y, XI)\n", + " -- Function File: YI = interp1 (..., METHOD)\n", + " -- Function File: YI = interp1 (..., EXTRAP)\n", + " -- Function File: YI = interp1 (..., \"left\")\n", + " -- Function File: YI = interp1 (..., \"right\")\n", + " -- Function File: PP = interp1 (..., \"pp\")\n", + "\n", + " One-dimensional interpolation.\n", + "\n", + " Interpolate input data to determine the value of YI at the points\n", + " XI. If not specified, X is taken to be the indices of Y ('1:length\n", + " (Y)'). If Y is a matrix or an N-dimensional array, the\n", + " interpolation is performed on each column of Y.\n", + "\n", + " The interpolation METHOD is one of:\n", + "\n", + " \"nearest\"\n", + " Return the nearest neighbor.\n", + "\n", + " \"previous\"\n", + " Return the previous neighbor.\n", + "\n", + " \"next\"\n", + " Return the next neighbor.\n", + "\n", + " \"linear\" (default)\n", + " Linear interpolation from nearest neighbors.\n", + "\n", + " \"pchip\"\n", + " Piecewise cubic Hermite interpolating\n", + " polynomial--shape-preserving interpolation with smooth first\n", + " derivative.\n", + "\n", + " \"cubic\"\n", + " Cubic interpolation (same as \"pchip\").\n", + "\n", + " \"spline\"\n", + " Cubic spline interpolation--smooth first and second\n", + " derivatives throughout the curve.\n", + "\n", + " Adding '*' to the start of any method above forces 'interp1' to\n", + " assume that X is uniformly spaced, and only 'X(1)' and 'X(2)' are\n", + " referenced. This is usually faster, and is never slower. The\n", + " default method is \"linear\".\n", + "\n", + " If EXTRAP is the string \"extrap\", then extrapolate values beyond\n", + " the endpoints using the current METHOD. If EXTRAP is a number,\n", + " then replace values beyond the endpoints with that number. When\n", + " unspecified, EXTRAP defaults to 'NA'.\n", + "\n", + " If the string argument \"pp\" is specified, then XI should not be\n", + " supplied and 'interp1' returns a piecewise polynomial object. This\n", + " object can later be used with 'ppval' to evaluate the\n", + " interpolation. There is an equivalence, such that 'ppval (interp1\n", + " (X, Y, METHOD, \"pp\"), XI) == interp1 (X, Y, XI, METHOD, \"extrap\")'.\n", + "\n", + " Duplicate points in X specify a discontinuous interpolant. There\n", + " may be at most 2 consecutive points with the same value. If X is\n", + " increasing, the default discontinuous interpolant is\n", + " right-continuous. If X is decreasing, the default discontinuous\n", + " interpolant is left-continuous. The continuity condition of the\n", + " interpolant may be specified by using the options \"left\" or \"right\"\n", + " to select a left-continuous or right-continuous interpolant,\n", + " respectively. Discontinuous interpolation is only allowed for\n", + " \"nearest\" and \"linear\" methods; in all other cases, the X-values\n", + " must be unique.\n", + "\n", + " An example of the use of 'interp1' is\n", + "\n", + " xf = [0:0.05:10];\n", + " yf = sin (2*pi*xf/5);\n", + " xp = [0:10];\n", + " yp = sin (2*pi*xp/5);\n", + " lin = interp1 (xp, yp, xf);\n", + " near = interp1 (xp, yp, xf, \"nearest\");\n", + " pch = interp1 (xp, yp, xf, \"pchip\");\n", + " spl = interp1 (xp, yp, xf, \"spline\");\n", + " plot (xf,yf,\"r\", xf,near,\"g\", xf,lin,\"b\", xf,pch,\"c\", xf,spl,\"m\",\n", + " xp,yp,\"r*\");\n", + " legend (\"original\", \"nearest\", \"linear\", \"pchip\", \"spline\");\n", + "\n", + " See also: pchip, spline, interpft, interp2, interp3, interpn.\n", + "\n", + "Additional help for built-in functions and operators is\n", + "available in the online version of the manual. Use the command\n", + "'doc ' to search the manual index.\n", + "\n", + "Help and information about Octave is also available on the WWW\n", + "at http://www.octave.org and via the help@octave.org\n", + "mailing list.\n" + ] + } + ], + "source": [ + "help interp1" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "'interp2' is a function from the file /usr/share/octave/4.0.0/m/general/interp2.m\n", + "\n", + " -- Function File: ZI = interp2 (X, Y, Z, XI, YI)\n", + " -- Function File: ZI = interp2 (Z, XI, YI)\n", + " -- Function File: ZI = interp2 (Z, N)\n", + " -- Function File: ZI = interp2 (Z)\n", + " -- Function File: ZI = interp2 (..., METHOD)\n", + " -- Function File: ZI = interp2 (..., METHOD, EXTRAP)\n", + "\n", + " Two-dimensional interpolation.\n", + "\n", + " Interpolate reference data X, Y, Z to determine ZI at the\n", + " coordinates XI, YI. The reference data X, Y can be matrices, as\n", + " returned by 'meshgrid', in which case the sizes of X, Y, and Z must\n", + " be equal. If X, Y are vectors describing a grid then 'length (X)\n", + " == columns (Z)' and 'length (Y) == rows (Z)'. In either case the\n", + " input data must be strictly monotonic.\n", + "\n", + " If called without X, Y, and just a single reference data matrix Z,\n", + " the 2-D region 'X = 1:columns (Z), Y = 1:rows (Z)' is assumed.\n", + " This saves memory if the grid is regular and the distance between\n", + " points is not important.\n", + "\n", + " If called with a single reference data matrix Z and a refinement\n", + " value N, then perform interpolation over a grid where each original\n", + " interval has been recursively subdivided N times. This results in\n", + " '2^N-1' additional points for every interval in the original grid.\n", + " If N is omitted a value of 1 is used. As an example, the interval\n", + " [0,1] with 'N==2' results in a refined interval with points at [0,\n", + " 1/4, 1/2, 3/4, 1].\n", + "\n", + " The interpolation METHOD is one of:\n", + "\n", + " \"nearest\"\n", + " Return the nearest neighbor.\n", + "\n", + " \"linear\" (default)\n", + " Linear interpolation from nearest neighbors.\n", + "\n", + " \"pchip\"\n", + " Piecewise cubic Hermite interpolating\n", + " polynomial--shape-preserving interpolation with smooth first\n", + " derivative.\n", + "\n", + " \"cubic\"\n", + " Cubic interpolation (same as \"pchip\").\n", + "\n", + " \"spline\"\n", + " Cubic spline interpolation--smooth first and second\n", + " derivatives throughout the curve.\n", + "\n", + " EXTRAP is a scalar number. It replaces values beyond the endpoints\n", + " with EXTRAP. Note that if EXTRAPVAL is used, METHOD must be\n", + " specified as well. If EXTRAP is omitted and the METHOD is\n", + " \"spline\", then the extrapolated values of the \"spline\" are used.\n", + " Otherwise the default EXTRAP value for any other METHOD is \"NA\".\n", + "\n", + " See also: interp1, interp3, interpn, meshgrid.\n", + "\n", + "Additional help for built-in functions and operators is\n", + "available in the online version of the manual. Use the command\n", + "'doc ' to search the manual index.\n", + "\n", + "Help and information about Octave is also available on the WWW\n", + "at http://www.octave.org and via the help@octave.org\n", + "mailing list.\n" + ] + } + ], + "source": [ + "help interp2" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t-1.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-0.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-3\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t3\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\tdata\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\tdata\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\n", + "\t\n", + "\tlinear\n", + "\n", + "\t\n", + "\t\tlinear\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tcubic spline\n", + "\n", + "\t\n", + "\t\tcubic spline\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tpiecewise cubic\n", + "\n", + "\t\n", + "\t\tpiecewise cubic\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "x=linspace(-pi,pi,9);\n", + "xi=linspace(-pi,pi,100);\n", + "y=sin(x);\n", + "yi_lin=interp1(x,y,xi,'linear');\n", + "yi_spline=interp1(x,y,xi,'spline'); \n", + "yi_cubic=interp1(x,y,xi,'cubic');\n", + "plot(x,y,'o',xi,yi_lin,xi,yi_spline,xi,yi_cubic)\n", + "axis([-pi,pi,-1.5,1.5])\n", + "legend('data','linear','cubic spline','piecewise cubic','Location','NorthWest')\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Example: Accelerate then hold velocity\n", + "\n", + "Here the time is given as vector t in seconds and the velocity is in mph. " + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t20\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t40\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t60\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t80\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t100\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t120\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t140\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t20\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t40\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t60\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t80\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t100\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t120\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\t\n", + "\t\tv (mph)\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\tt (s)\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\tdata\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\tdata\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\n", + "\t\n", + "\tlinear\n", + "\n", + "\t\n", + "\t\tlinear\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tcubic spline\n", + "\n", + "\t\n", + "\t\tcubic spline\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tpiecewise cubic\n", + "\n", + "\t\n", + "\t\tpiecewise cubic\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "t=[0 20 40 56 68 80 84 96 104 110]';\n", + "v=[0 20 20 38 80 80 100 100 125 125]';\n", + "tt=linspace(0,110)';\n", + "v_lin=interp1(t,v,tt);\n", + "v_spl=interp1(t,v,tt,'spline');\n", + "v_cub=interp1(t,v,tt,'cubic');\n", + "\n", + "plot(t,v,'o',tt,v_lin,tt,v_spl,tt,v_cub)\n", + "xlabel('t (s)')\n", + "ylabel('v (mph)')\n", + "legend('data','linear','cubic spline','piecewise cubic','Location','NorthWest')" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t-4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t-2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t6\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t8\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t20\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t40\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t60\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t80\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t100\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t120\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\t\n", + "\t\tdv/dt (mph/s)\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\tt (s)\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\tlinear\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\tlinear\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tcubic spline\n", + "\n", + "\t\n", + "\t\tcubic spline\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tpiecewise cubic\n", + "\n", + "\t\n", + "\t\tpiecewise cubic\n", + "\t\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "t=[0 20 40 56 68 80 84 96 104 110]';\n", + "v=[0 20 20 38 80 80 100 100 125 125]';\n", + "tt=linspace(0,110)';\n", + "v_lin=interp1(t,v,tt);\n", + "v_spl=interp1(t,v,tt,'spline');\n", + "v_cub=interp1(t,v,tt,'cubic');\n", + "\n", + "\n", + "plot(tt(2:end),diff(v_lin)./diff(tt),tt(2:end),diff(v_spl)./diff(tt),tt(2:end),diff(v_cub)./diff(tt))\n", + "xlabel('t (s)')\n", + "ylabel('dv/dt (mph/s)')\n", + "legend('linear','cubic spline','piecewise cubic','Location','NorthWest')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Choose spline wisely\n", + "\n", + "For example of sin(x), not very important\n", + "\n", + "For stop-and-hold examples, the $C^{2}$-continuity should not be preserved. You don't need smooth curves.\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Numerical Integration\n", + "\n", + "A definite integral is defined by \n", + "\n", + "$I=\\int_{a}^{b}f(x)dx$\n", + "\n", + "To determine the mass of an object with varying density, you can perform a summation\n", + "\n", + "mass=$\\sum_{i=1}^{n}\\rho_{i}\\Delta V_{i}$\n", + "\n", + "or taking the limit as $\\Delta V \\rightarrow dV=dxdydz$\n", + "\n", + "mass=$\\int_{0}^{h}\\int_{0}^{w}\\int_{0}^{l}\\rho(x,y,z)dxdydz$\n", + "\n", + "## Newton-Cotes Formulas\n", + "\n", + "$I=\\int_{a}^{b}f(x)dx=\\int_{a}^{b}f_{n}(x)dx$\n", + "\n", + "where $f_{n}$ is an n$^{th}$-order polynomial approximation of f(x)\n", + "\n", + "## First-Order: Trapezoidal Rule\n", + "\n", + "$I=\\int_{a}^{b}f(x)dx\\approx \\int_{a}^{b}\\left(f(a)+\\frac{f(b)-f(a)}{b-a}(x-a)\\right)dx$\n", + "\n", + "$I\\approx(b-a)\\frac{f(a)+f(b)}{2}$" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "I_trap = 0.78540\n", + "I_act = 1.00000\n" + ] + }, + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.6\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.8\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t2.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t3\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t3.5\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\n", + "\tgnuplot_plot_1a\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tgnuplot_plot_2a\n", + "\n", + "\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "x=linspace(0,pi)';\n", + "plot(x,sin(x),[0,pi/2],sin([0,pi/2]))\n", + "I_trap=mean(sin(([0,pi/2]))*(diff([0,pi/2])))\n", + "I_act = -(cos(pi/2)-cos(0))\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Improve estimate with more points\n", + "\n", + "$I=\\int_{a}^{b}f(x)dx=\\int_{a}^{a+\\Delta x}f(x)dx+\\int_{a+\\Delta x}^{a+2\\Delta x}f(x)dx+ \\cdots \\int_{b-\\Delta x}^{b}f(x)dx$\n", + "\n", + "$I\\approx\\Delta x\\frac{f(a)+f(a+\\Delta x)}{2}+\\Delta x\\frac{f(a+\\Delta x)+f(a+2\\Delta x)}{2}\n", + "+\\cdots \\Delta x\\frac{f(b-\\Delta x)+f(b)}{2}$\n", + "\n", + "$I\\approx \\frac{\\Delta x}{2}\\left(f(a)+2\\sum_{i=1}^{n-1}f(a+i\\Delta x) +f(b)\\right)$" + ] + }, + { + "cell_type": "code", + "execution_count": 48, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "For 5 steps\n", + "trapezoid approximation of integral is 0.79 \n", + " actual integral is 1.00\n" + ] + } + ], + "source": [ + "N=5;\n", + "I_trap=trap(@(x) sin(x),0,pi/2,N);\n", + "fprintf('For %i steps\\ntrapezoid approximation of integral is %1.2f \\n actual integral is %1.2f',N,I_trap,I_act)" + ] + }, + { + "cell_type": "code", + "execution_count": 43, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "Gnuplot\n", + "Produced by GNUPLOT 5.0 patchlevel 3 \n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\t\n", + "\t \n", + "\t \n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\t\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.4\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.6\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.8\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t0.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t1.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t2\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t2.5\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t3\n", + "\t\n", + "\n", + "\n", + "\t\t\n", + "\t\t3.5\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\t\n", + "\n", + "\n", + "\tgnuplot_plot_1a\n", + "\n", + "\n", + "\n", + "\t\n", + "\t\n", + "\tgnuplot_plot_2a\n", + "\n", + "\t \n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\t\n", + "\n", + "\t\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "\n", + "plot(x,sin(x),linspace(0,pi/2,N),sin(linspace(0,pi/2,N)),'o')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Increase accuracy = Increase polynomial order\n", + "\n", + "### Simpson's Rules\n", + "\n", + "When integrating f(x) and using a second order polynomial, this is known as **Simpson's 1/3 Rule**\n", + "\n", + "$I=\\frac{h}{3}(f(x_{0})+4f(x_{1})+f(x_{2}))$\n", + "\n", + "where a=$x_{0}$, b=$x_{2}$, and $x_{1}=\\frac{a+b}{2}$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This can be used with n=3 or multiples of 2 intervals\n", + "\n", + "$I=\\int_{x_{0}}^{x_{2}}f(x)dx+\\int_{x_{2}}^{x_{4}}f(x)dx+\\cdots +\\int_{x_{n-2}}^{x_{n}}f(x)dx$\n", + "\n", + "$I=(b-a)\\frac{f(x_{0})+4\\sum_{i=1,3,5}^{n-1}f(x_{i})+2\\sum_{i=2,4,6}^{n-2}f(x_{i})+f(x_{n})}{3n}$" + ] + }, + { + "cell_type": "code", + "execution_count": 68, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ans = 1.6235\n", + "Is_1_3 = 1.0023\n" + ] + } + ], + "source": [ + "f=@(x) 0.2+25*x-200*x.^2+675*x.^3-900*x.^4+400*x.^5;\n", + "simpson3(f,0,0.8,4)\n", + "Is_1_3=simpson3(@(x) sin(x),0,pi/2,2)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## General Newton-Cotes formulae\n", + "\n", + "![Newton-Cotes Table](newton_cotes.png)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Octave", + "language": "octave", + "name": "octave" + }, + "language_info": { + "file_extension": ".m", + "help_links": [ + { + "text": "MetaKernel Magics", + "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md" + } + ], + "mimetype": "text/x-octave", + "name": "octave", + "version": "0.19.14" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/lecture_19/newton_cotes.png b/lecture_19/newton_cotes.png new file mode 100644 index 0000000..b734d11 Binary files /dev/null and b/lecture_19/newton_cotes.png differ diff --git a/lecture_19/q1.png b/lecture_19/q1.png new file mode 100644 index 0000000..44ec9da Binary files /dev/null and b/lecture_19/q1.png differ diff --git a/lecture_19/q2.png b/lecture_19/q2.png new file mode 100644 index 0000000..fb27a54 Binary files /dev/null and b/lecture_19/q2.png differ diff --git a/lecture_19/q3.png b/lecture_19/q3.png new file mode 100644 index 0000000..62c4402 Binary files /dev/null and b/lecture_19/q3.png differ diff --git a/lecture_19/q4.png b/lecture_19/q4.png new file mode 100644 index 0000000..5457483 Binary files /dev/null and b/lecture_19/q4.png differ diff --git a/lecture_19/simpson3.m b/lecture_19/simpson3.m new file mode 100644 index 0000000..2ae0c81 --- /dev/null +++ b/lecture_19/simpson3.m @@ -0,0 +1,23 @@ +function I = simpson3(func,a,b,n,varargin) +% simpson3: composite simpson's 1/3 rule +% I = simpson3(func,a,b,n,pl,p2,...): +% composite trapezoidal rule +% input: +% func = name of function to be integrated +% a, b = integration limits +% n = number of segments (default = 100) +% pl,p2,... = additional parameters used by func +% output: +% I = integral estimate +if nargin<3,error('at least 3 input arguments required'),end +if ~(b>a),error('upper bound must be greater than lower'),end +if nargin<4|isempty(n),n=100;end +x = a; h = (b - a)/n; + +xvals=linspace(a,b,n+1); +fvals=func(xvals,varargin{:}); +s=fvals(1); +s = s + 4*sum(fvals(2:2:end-1)); +s = s + 2*sum(fvals(3:2:end-2)); +s = s + fvals(end); +I = (b - a) * s/(3*n); diff --git a/lecture_19/trap.m b/lecture_19/trap.m new file mode 100644 index 0000000..85b8685 --- /dev/null +++ b/lecture_19/trap.m @@ -0,0 +1,22 @@ +function I = trap(func,a,b,n,varargin) +% trap: composite trapezoidal rule quadrature +% I = trap(func,a,b,n,pl,p2,...): +% composite trapezoidal rule +% input: +% func = name of function to be integrated +% a, b = integration limits +% n = number of segments (default = 100) +% pl,p2,... = additional parameters used by func +% output: +% I = integral estimate +if nargin<3,error('at least 3 input arguments required'),end +if ~(b>a),error('upper bound must be greater than lower'),end +if nargin<4|isempty(n),n=100;end + +x = a; h = (b - a)/n; +xvals=linspace(a,b,n); +fvals=func(xvals,varargin{:}); +s=func(a,varargin{:}); +s = s + 2*sum(fvals(2:n-1)); +s = s + func(b,varargin{:}); +I = (b - a) * s/(2*n);