Interpolation and ApproximationMirko Navara (navara@math.feld.cvut.cz) & Ale\305\241 N\304\233me\304\215ek (nemecek@math.feld.cvut.cz)Program descriptionThis worksheet implementsPolynomial interpolation - Lagrange's classical formulaPolynomial interpolation - Newtons's 1st formulaPolynomial interpolation - Newtons's 2nd formulaPolynomial interpolation - standard function interpSpline interpolationPlain interpolationLeast Squares Approximation (LSM)Taylor's approximation - standard function taylorGlobal variables:x, y ... lists of coordinates ([x[i], y[i]] - point) n ... number of pointsf ... given functionphi ... interpolating/approximating functionpsi[i] $i=1..k ... auxiliary functionsc[i] $i=1..k ... coefficientsa..b ... range of interpolation/aproximation (or plot interval)eps ... given accuracy
Do NOT modify commands marked by ## !Usually do not modify commands in standard math input format.Each section ends with links to other parts of the program.Program initializationRun this part only once (possibly after LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYoLUkjbWlHRiQ2JlEocmVzdGFydEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUrZXhlY3V0YWJsZUdGMS8lLG1hdGh2YXJpYW50R1EnaXRhbGljRicvJSVzaXplR1EjMTJGJy8lK2ZvcmVncm91bmRHUShbMCwwLDBdRicvJStiYWNrZ3JvdW5kR1EuWzI1NSwyNTUsMjU1XUYnRjIvRjVRJ25vcm1hbEYn).with(plots, display, setoptions): ##with(Student[LinearAlgebra], [LeastSquaresPlot, LeastSquares]): ##with(linalg, leastsqrs):with(Optimization, Maximize, Minimize): ##Procedure calculating coefficients:koeficienty:= proc() local solver, substituce; ## unassign(''c[i]'' $i=1..k): # unassign names if (metoda="MNC" and MNCoption="new")
then solver:='LinearAlgebra[LeastSquares]'
elif (metoda="MNC" and MNCoption="old") then solver:='linal[leastsqrs]'
else solver:='fsolve'
end if; # use fsolve/solve substituce:=solver({'phi(x[i])=y[i]' $i=1..n},
{'c[i]' $i=1..k}); assign(substituce); substituce end: ##plots[setoptions](thickness=1, symbol=solidcircle, symbolsize=12, titlefont=[TIMES, ROMAN, 15], legendstyle=[font=[TIMES, ROMAN, 15]]);AssignmentChoice of a) interpolation/approximation of a given function b) fitting a set of pointsa) Given function and interval (you choose points)Option 1: Read data from fileIf the file is not in the working directory, specify the path and the name of the file (e.g. "~/nm/myjob/uloha101E.txt").read("uloha1??E.txt");eps:=epsilon:Option 2: Manual inputIf you change the assignment, you must define f as a function/procedure, e.g. f := t -> 3 + sin(t), not an expression.f:=sin;a:=0:b:=5:Common part:LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYqLUkjbWlHRiQ2I1EhRictRiM2Ky1GLDYmUSlpbnRlcnZhbEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUrZXhlY3V0YWJsZUdGNi8lLG1hdGh2YXJpYW50R1EnaXRhbGljRictSSNtb0dGJDYuUSM6PUYnRjcvRjpRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0ZELyUpc3RyZXRjaHlHRkQvJSpzeW1tZXRyaWNHRkQvJShsYXJnZW9wR0ZELyUubW92YWJsZWxpbWl0c0dGRC8lJ2FjY2VudEdGRC8lJ2xzcGFjZUdRLDAuMjc3Nzc3OGVtRicvJSdyc3BhY2VHRlMtRiM2Ki1GLDYmUSJhRidGNEY3RjktRj02LlEjLi5GJ0Y3RkBGQkZFRkdGSUZLRk1GTy9GUlEsMC4yMjIyMjIyZW1GJy9GVVEmMC4wZW1GJy1GLDYmUSJiRidGNEY3RjkvJSVzaXplR1EjMTJGJy8lK2ZvcmVncm91bmRHUShbMCwwLDBdRicvJStiYWNrZ3JvdW5kR1EuWzI1NSwyNTUsMjU1XUYnRjdGQEYrRl9vRmJvRmVvRjdGQEYrRl9vRmJvRmVvRjdGQA==Number of points:n:=5;JSFHOption 1: Equidistant sequence of points Coordinates x (real values induce floating point arithmetic).LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYqLUkjbWlHRiQ2I1EhRictRiM2Ki1GLDYmUSJ4RicvJSdpdGFsaWNHUSV0cnVlRicvJStleGVjdXRhYmxlR0Y2LyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi5RIzo9RidGNy9GOlEnbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRkQvJSlzdHJldGNoeUdGRC8lKnN5bW1ldHJpY0dGRC8lKGxhcmdlb3BHRkQvJS5tb3ZhYmxlbGltaXRzR0ZELyUnYWNjZW50R0ZELyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGUy1JKG1mZW5jZWRHRiQ2Jy1GIzYqRistRiM2K0YrLUYjNictRiw2JlEiYUYnRjRGN0Y5LUY9Ni5RIitGJ0Y3RkBGQkZFRkdGSUZLRk1GTy9GUlEsMC4yMjIyMjIyZW1GJy9GVUZgby1JJm1mcmFjR0YkNigtRiM2Ki1GVzYlLUYjNictRiw2JlEiYkYnRjRGN0Y5LUY9Ni5RKCZtaW51cztGJ0Y3RkBGQkZFRkdGSUZLRk1GT0Zfb0Zhb0ZpbkY3RkBGN0ZALUY9Ni5RMSZJbnZpc2libGVUaW1lcztGJ0Y3RkBGQkZFRkdGSUZLRk1GTy9GUlEmMC4wZW1GJy9GVUZlcC1GVzYlLUYjNictRiw2JlEiaUYnRjRGN0Y5Rl5wLUkjbW5HRiQ2JVEiMUYnRjdGQEY3RkBGN0ZALyUlc2l6ZUdRIzEyRicvJStmb3JlZ3JvdW5kR1EoWzAsMCwwXUYnLyUrYmFja2dyb3VuZEdRLlsyNTUsMjU1LDI1NV1GJ0Y3RkAtRiM2KkYrLUYjNictRiw2JlEibkYnRjRGN0Y5Rl5wRl5xRjdGQEYrRmJxRmVxRmhxRjdGQC8lLmxpbmV0aGlja25lc3NHRmFxLyUrZGVub21hbGlnbkdRJ2NlbnRlckYnLyUpbnVtYWxpZ25HRmZyLyUpYmV2ZWxsZWRHRkRGN0ZALUY9Ni5RIiRGJ0Y3RkBGQkZFRkdGSUZLRk1GT0ZkcEZmcC1GVzYlLUYjNitGW3EtRj02LlEiPUYnRjdGQEZCRkVGR0ZJRktGTUZPRlFGVC1GIzYqRl5xLUY9Ni5RIy4uRidGN0ZARkJGRUZHRklGS0ZNRk9GX29GZnBGX3JGYnFGZXFGaHFGN0ZARitGYnFGZXFGaHFGN0ZARjdGQEZicUZlcUZocUY3RkBGK0ZicUZlcUZocUY3RkBGN0ZALyUlb3BlbkdRIltGJy8lJmNsb3NlR1EiXUYnRmJxRmVxRmhxRjdGQEYrRmJxRmVxRmhxRjdGQA==Coordinates y:LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYqLUkjbWlHRiQ2I1EhRictRiM2Ky1GLDYmUSJ5RicvJSdpdGFsaWNHUSV0cnVlRicvJStleGVjdXRhYmxlR0Y2LyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi5RIzo9RidGNy9GOlEnbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRkQvJSlzdHJldGNoeUdGRC8lKnN5bW1ldHJpY0dGRC8lKGxhcmdlb3BHRkQvJS5tb3ZhYmxlbGltaXRzR0ZELyUnYWNjZW50R0ZELyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGUy1GIzYqLUYsNiZRJG1hcEYnRjRGN0Y5LUY9Ni5RMCZBcHBseUZ1bmN0aW9uO0YnRjdGQEZCRkVGR0ZJRktGTUZPL0ZSUSYwLjBlbUYnL0ZVRmluLUkobWZlbmNlZEdGJDYlLUYjNiotRiw2JlEiZkYnRjRGN0Y5LUY9Ni5RIixGJ0Y3RkBGQi9GRkY2RkdGSUZLRk1GT0Zobi9GVVEsMC4zMzMzMzMzZW1GJy1GLDYmUSJ4RidGNEY3RjkvJSVzaXplR1EjMTJGJy8lK2ZvcmVncm91bmRHUShbMCwwLDBdRicvJStiYWNrZ3JvdW5kR1EuWzI1NSwyNTUsMjU1XUYnRjdGQEY3RkBGXHBGX3BGYnBGN0ZARitGXHBGX3BGYnBGN0ZARitGXHBGX3BGYnBGN0ZAJSFHOption 2: Cosine sequence of pointsLUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYnLUkjbWlHRiQ2I1EhRictRiM2Jy1GLDYmUSJ4RicvJSdpdGFsaWNHUSV0cnVlRicvJStleGVjdXRhYmxlR0Y2LyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi5RIzo9RidGNy9GOlEnbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRkQvJSlzdHJldGNoeUdGRC8lKnN5bW1ldHJpY0dGRC8lKGxhcmdlb3BHRkQvJS5tb3ZhYmxlbGltaXRzR0ZELyUnYWNjZW50R0ZELyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGUy1JKG1mZW5jZWRHRiQ2Jy1GIzYnRistRiM2KEYrLUYjNictSSZtZnJhY0dGJDYoLUYjNidGKy1GIzYnLUYsNiZRImFGJ0Y0RjdGOS1GPTYuUSIrRidGN0ZARkJGRUZHRklGS0ZNRk8vRlJRLDAuMjIyMjIyMmVtRicvRlVGZ28tRiw2JlEiYkYnRjRGN0Y5RjdGQEYrRjdGQC1GIzYlLUkjbW5HRiQ2JVEiMkYnRjdGQEY3RkAvJS5saW5ldGhpY2tuZXNzR1EiMUYnLyUrZGVub21hbGlnbkdRJ2NlbnRlckYnLyUpbnVtYWxpZ25HRmdwLyUpYmV2ZWxsZWRHRkRGY28tRmpuNigtRiM2KC1GVzYlLUYjNidGYG8tRj02LlEoJm1pbnVzO0YnRjdGQEZCRkVGR0ZJRktGTUZPRmZvRmhvRmlvRjdGQEY3RkAtRj02LlExJkludmlzaWJsZVRpbWVzO0YnRjdGQEZCRkVGR0ZJRktGTUZPL0ZSUSYwLjBlbUYnL0ZVRltyLUYjNictRiw2JlEkY29zRicvRjVGREY3RkAtRj02LlEwJkFwcGx5RnVuY3Rpb247RidGN0ZARkJGRUZHRklGS0ZNRk9GanFGXHItRlc2JS1GIzYlLUZqbjYoLUYjNictRlc2JS1GIzYnLUYsNiZRImlGJ0Y0RjdGOUZkcS1Gam42KC1GIzYlLUZfcDYlRmRwRjdGQEY3RkBGXHBGYnBGZXBGaHBGanBGN0ZARjdGQEZncS1GLDYmUSUmcGk7RidGYnJGN0ZARjdGQC1GIzYlLUYsNiZRIm5GJ0Y0RjdGOUY3RkBGYnBGZXBGaHBGanBGN0ZARjdGQEY3RkBGK0Y3RkBGXHBGYnBGZXBGaHBGanBGN0ZALUY9Ni5RIiRGJ0Y3RkBGQkZFRkdGSUZLRk1GT0ZqcUZcci1GVzYlLUYjNihGYnMtRj02LlEiPUYnRjdGQEZCRkVGR0ZJRktGTUZPRlFGVC1GIzYnRmlzLUY9Ni5RIy4uRidGN0ZARkJGRUZHRklGS0ZNRk9GZm9GXHJGYHRGN0ZARitGN0ZARjdGQEY3RkBGK0Y3RkBGN0ZALyUlb3BlbkdRIltGJy8lJmNsb3NlR1EiXUYnRjdGQEYrRjdGQA==Coordinates y:LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYnLUkjbWlHRiQ2I1EhRictRiM2KC1GLDYmUSJ5RicvJSdpdGFsaWNHUSV0cnVlRicvJStleGVjdXRhYmxlR0Y2LyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi5RIzo9RidGNy9GOlEnbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRkQvJSlzdHJldGNoeUdGRC8lKnN5bW1ldHJpY0dGRC8lKGxhcmdlb3BHRkQvJS5tb3ZhYmxlbGltaXRzR0ZELyUnYWNjZW50R0ZELyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGUy1GIzYnLUYsNiZRJG1hcEYnRjRGN0Y5LUY9Ni5RMCZBcHBseUZ1bmN0aW9uO0YnRjdGQEZCRkVGR0ZJRktGTUZPL0ZSUSYwLjBlbUYnL0ZVRmluLUkobWZlbmNlZEdGJDYlLUYjNictRiw2JlEiZkYnRjRGN0Y5LUY9Ni5RIixGJ0Y3RkBGQi9GRkY2RkdGSUZLRk1GT0Zobi9GVVEsMC4zMzMzMzMzZW1GJy1GLDYmUSJ4RidGNEY3RjlGN0ZARjdGQEY3RkBGK0Y3RkBGK0Y3RkA=Option 3: Customized sequence of pointsx:=[1.+i/3 $i=0..5, 3, 4, 4+i/2 $i=1..4];LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYnLUkjbWlHRiQ2I1EhRictRiM2KC1GLDYmUSJ5RicvJSdpdGFsaWNHUSV0cnVlRicvJStleGVjdXRhYmxlR0Y2LyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi5RIzo9RidGNy9GOlEnbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRkQvJSlzdHJldGNoeUdGRC8lKnN5bW1ldHJpY0dGRC8lKGxhcmdlb3BHRkQvJS5tb3ZhYmxlbGltaXRzR0ZELyUnYWNjZW50R0ZELyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGUy1GIzYnLUYsNiZRJG1hcEYnRjRGN0Y5LUY9Ni5RMCZBcHBseUZ1bmN0aW9uO0YnRjdGQEZCRkVGR0ZJRktGTUZPL0ZSUSYwLjBlbUYnL0ZVRmluLUkobWZlbmNlZEdGJDYlLUYjNictRiw2JlEiZkYnRjRGN0Y5LUY9Ni5RIixGJ0Y3RkBGQi9GRkY2RkdGSUZLRk1GT0Zobi9GVVEsMC4zMzMzMzMzZW1GJy1GLDYmUSJ4RidGNEY3RjlGN0ZARjdGQEY3RkBGK0Y3RkBGK0Y3RkA=LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYnLUkjbWlHRiQ2I1EhRictRiM2KC1GLDYmUSJuRicvJSdpdGFsaWNHUSV0cnVlRicvJStleGVjdXRhYmxlR0Y2LyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi5RIzo9RidGNy9GOlEnbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRkQvJSlzdHJldGNoeUdGRC8lKnN5bW1ldHJpY0dGRC8lKGxhcmdlb3BHRkQvJS5tb3ZhYmxlbGltaXRzR0ZELyUnYWNjZW50R0ZELyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGUy1GIzYnLUYsNiZRJW5vcHNGJ0Y0RjdGOS1GPTYuUTAmQXBwbHlGdW5jdGlvbjtGJ0Y3RkBGQkZFRkdGSUZLRk1GTy9GUlEmMC4wZW1GJy9GVUZpbi1JKG1mZW5jZWRHRiQ2JS1GIzYlLUYsNiZRInhGJ0Y0RjdGOUY3RkBGN0ZARjdGQEYrRjdGQEYrRjdGQA==Common option: Convert symbolic values to floating points numbersLUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYnLUkjbWlHRiQ2I1EhRictRiM2KC1GLDYmUSJ4RicvJSdpdGFsaWNHUSV0cnVlRicvJStleGVjdXRhYmxlR0Y2LyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi5RIzo9RidGNy9GOlEnbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRkQvJSlzdHJldGNoeUdGRC8lKnN5bW1ldHJpY0dGRC8lKGxhcmdlb3BHRkQvJS5tb3ZhYmxlbGltaXRzR0ZELyUnYWNjZW50R0ZELyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGUy1GIzYnLUYsNiZRJmV2YWxmRidGNEY3RjktRj02LlEwJkFwcGx5RnVuY3Rpb247RidGN0ZARkJGRUZHRklGS0ZNRk8vRlJRJjAuMGVtRicvRlVGaW4tSShtZmVuY2VkR0YkNiUtRiM2JUYxRjdGQEY3RkBGN0ZARitGN0ZARitGN0ZALUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYnLUkjbWlHRiQ2I1EhRictRiM2KC1GLDYmUSJ5RicvJSdpdGFsaWNHUSV0cnVlRicvJStleGVjdXRhYmxlR0Y2LyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi5RIzo9RidGNy9GOlEnbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRkQvJSlzdHJldGNoeUdGRC8lKnN5bW1ldHJpY0dGRC8lKGxhcmdlb3BHRkQvJS5tb3ZhYmxlbGltaXRzR0ZELyUnYWNjZW50R0ZELyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGUy1GIzYnLUYsNiZRJmV2YWxmRidGNEY3RjktRj02LlEwJkFwcGx5RnVuY3Rpb247RidGN0ZARkJGRUZHRklGS0ZNRk8vRlJRJjAuMGVtRicvRlVGaW4tSShtZmVuY2VkR0YkNiUtRiM2JUYxRjdGQEY3RkBGN0ZARitGN0ZARitGN0ZAJSFHGo to function definition, table of points, choice of approximation: Lagrange, Newton 1, Newton 2, interp, spline, plain, , Taylor,different forms of output, graphical outputs for a given function , table of points.b) Given table of points onlyOption 1: Read data from fileIf the file is not in the working directory, specify the path and the name of the file (e.g. "~/nm/myjob/uloha101E.txt").read("uloha2??E.txt");Option 2: Manual inputx:=[$1.0..4.0,6.0]; y:=[1.,2.,1.,1.,1.];Common part:LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYqLUkjbWlHRiQ2I1EhRictRiM2Ky1GLDYmUSJuRicvJSdpdGFsaWNHUSV0cnVlRicvJStleGVjdXRhYmxlR0Y2LyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi5RIzo9RidGNy9GOlEnbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRkQvJSlzdHJldGNoeUdGRC8lKnN5bW1ldHJpY0dGRC8lKGxhcmdlb3BHRkQvJS5tb3ZhYmxlbGltaXRzR0ZELyUnYWNjZW50R0ZELyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGUy1GIzYqLUYsNiZRJW5vcHNGJ0Y0RjdGOS1GPTYuUTAmQXBwbHlGdW5jdGlvbjtGJ0Y3RkBGQkZFRkdGSUZLRk1GTy9GUlEmMC4wZW1GJy9GVUZpbi1JKG1mZW5jZWRHRiQ2JS1GIzYoLUYsNiZRInhGJ0Y0RjdGOS8lJXNpemVHUSMxMkYnLyUrZm9yZWdyb3VuZEdRKFswLDAsMF1GJy8lK2JhY2tncm91bmRHUS5bMjU1LDI1NSwyNTVdRidGN0ZARjdGQEZjb0Zmb0Zpb0Y3RkBGK0Zjb0Zmb0Zpb0Y3RkBGK0Zjb0Zmb0Zpb0Y3RkA=LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYnLUkjbWlHRiQ2I1EhRictRiM2KC1GLDYmUSJhRicvJSdpdGFsaWNHUSV0cnVlRicvJStleGVjdXRhYmxlR0Y2LyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi5RIzo9RidGNy9GOlEnbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRkQvJSlzdHJldGNoeUdGRC8lKnN5bW1ldHJpY0dGRC8lKGxhcmdlb3BHRkQvJS5tb3ZhYmxlbGltaXRzR0ZELyUnYWNjZW50R0ZELyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGUy1GIzYnLUYsNiZRJG1pbkYnL0Y1RkRGN0ZALUY9Ni5RMCZBcHBseUZ1bmN0aW9uO0YnRjdGQEZCRkVGR0ZJRktGTUZPL0ZSUSYwLjBlbUYnL0ZVRmpuLUkobWZlbmNlZEdGJDYlLUYjNidGKy1GIzYnLUYsNiZRI29wRidGNEY3RjlGZm4tRl1vNiUtRiM2JS1GLDYmUSJ4RidGNEY3RjlGN0ZARjdGQEY3RkBGK0Y3RkBGN0ZARjdGQEYrRjdGQEYrRjdGQA==LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYnLUkjbWlHRiQ2I1EhRictRiM2KC1GLDYmUSJiRicvJSdpdGFsaWNHUSV0cnVlRicvJStleGVjdXRhYmxlR0Y2LyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1JI21vR0YkNi5RIzo9RidGNy9GOlEnbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRkQvJSlzdHJldGNoeUdGRC8lKnN5bW1ldHJpY0dGRC8lKGxhcmdlb3BHRkQvJS5tb3ZhYmxlbGltaXRzR0ZELyUnYWNjZW50R0ZELyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGUy1GIzYnLUYsNiZRJG1heEYnL0Y1RkRGN0ZALUY9Ni5RMCZBcHBseUZ1bmN0aW9uO0YnRjdGQEZCRkVGR0ZJRktGTUZPL0ZSUSYwLjBlbUYnL0ZVRmpuLUkobWZlbmNlZEdGJDYlLUYjNidGKy1GIzYnLUYsNiZRI29wRidGNEY3RjlGZm4tRl1vNiUtRiM2JS1GLDYmUSJ4RidGNEY3RjlGN0ZARjdGQEY3RkBGK0Y3RkBGN0ZARjdGQEYrRjdGQEYrRjdGQA==LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYnLUkjbWlHRiQ2I1EhRictRiM2KC1GLDYmUSlpbnRlcnZhbEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUrZXhlY3V0YWJsZUdGNi8lLG1hdGh2YXJpYW50R1EnaXRhbGljRictSSNtb0dGJDYuUSM6PUYnRjcvRjpRJ25vcm1hbEYnLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0ZELyUpc3RyZXRjaHlHRkQvJSpzeW1tZXRyaWNHRkQvJShsYXJnZW9wR0ZELyUubW92YWJsZWxpbWl0c0dGRC8lJ2FjY2VudEdGRC8lJ2xzcGFjZUdRLDAuMjc3Nzc3OGVtRicvJSdyc3BhY2VHRlMtRiM2Jy1GLDYmUSJhRidGNEY3RjktRj02LlEjLi5GJ0Y3RkBGQkZFRkdGSUZLRk1GTy9GUlEsMC4yMjIyMjIyZW1GJy9GVVEmMC4wZW1GJy1GLDYmUSJiRidGNEY3RjlGN0ZARitGN0ZARitGN0ZALUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYnLUkjbWlHRiQ2I1EhRictRiM2KUYrLUYjNi8tSSNtb0dGJDYwUSNpZkYnLyUlYm9sZEdRJXRydWVGJy8lK2V4ZWN1dGFibGVHRjkvJSxtYXRodmFyaWFudEdRJWJvbGRGJy8lK2ZvbnR3ZWlnaHRHRj4vJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRkMvJSlzdHJldGNoeUdGQy8lKnN5bW1ldHJpY0dGQy8lKGxhcmdlb3BHRkMvJS5tb3ZhYmxlbGltaXRzR0ZDLyUnYWNjZW50R0ZDLyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdGUi1JJ21zcGFjZUdGJDYmLyUnaGVpZ2h0R1EmMC4wZXhGJy8lJndpZHRoR1EmMC41ZW1GJy8lJmRlcHRoR0ZaLyUqbGluZWJyZWFrR1ElYXV0b0YnLUYjNigtRiw2JlEibkYnLyUnaXRhbGljR0Y5RjovRj1RJ2l0YWxpY0YnLUY0Ni5RIj1GJ0Y6L0Y9USdub3JtYWxGJ0ZBRkRGRkZIRkpGTEZOL0ZRUSwwLjI3Nzc3NzhlbUYnL0ZURlxwLUYjNictRiw2JlElbm9wc0YnRmJvRjpGZG8tRjQ2LlEwJkFwcGx5RnVuY3Rpb247RidGOkZpb0ZBRkRGRkZIRkpGTEZORlBGUy1JKG1mZW5jZWRHRiQ2JS1GIzYlLUYsNiZRInlGJ0Zib0Y6RmRvRjpGaW9GOkZpb0Y6RmlvRitGOkZpb0YrRlUtRjQ2MFEldGhlbkYnRjdGOkY8Rj9GQUZERkZGSEZKRkxGTkZQRlMtRlY2JkZYRmVuRmhuL0Zbb1E2aW5jcmVhc2VpbmRlbnRuZXdsaW5lRictRiM2Jy1GLDYmUSdscHJpbnRGJ0Zib0Y6RmRvRmNwLUZncDYlLUYjNilGKy1JI21zR0YkNiNRMk51bWJlcn5vZn5wb2ludHM6RictRjQ2LlEiLEYnRjpGaW9GQS9GRUY5RkZGSEZKRkxGTkZQL0ZUUSwwLjMzMzMzMzNlbUYnLUYjNihGKy1GIzYnLUY0Ni5RIidGJ0Y6RmlvRkFGREZGRkhGSkZMRk4vRlFRLDAuMTExMTExMWVtRidGU0Zfb0Zcc0Y6RmlvRmZvRl9vRjpGaW9GK0Y6RmlvRjpGaW9GOkZpb0YrLUZWNiZGWEZlbkZobi9GW29RNmRlY3JlYXNlaW5kZW50bmV3bGluZUYnRitGOkZpb0YrLUYjNiotRjQ2MFElZWxzZUYnRjdGOkY8Rj9GQUZERkZGSEZKRkxGTkZQRlNGYXEtRiM2J0ZncUZjcC1GZ3A2JS1GIzYnRistRl9yNiNRTEVycm9yOn5kaW1lbnNpb25zfm9mfmxpc3Rzfngsfnl+YXJlfnVuZXF1YWxGJ0YrRjpGaW9GOkZpb0Y6RmlvRitGYXNGK0Y6RmlvLUY0NjBRJ2VuZH5pZkYnRjdGOkY8Rj9GQUZERkZGSEZKRkxGTkZQRlNGOkZpb0YrRjpGaW8=Go to function definition, table of points, choice of approximation: Lagrange, Newton 1, Newton 2, interp, spline, plain, , Taylor,different forms of output, graphical outputs for a given function , table of points.Interpolation/Approximation methodsCHOICE OF METHODS OF INTERPOLATION/APPROXIMATION (or make your own):1) Polynomial interpolation - Lagrange's classical formulametoda:="LAGRANGE": k:=n: ##unassign('i','j'): ##LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUkjbWlHRiQ2I1EhRictRiM2JkYrLUYjNiYtRiw2JVEkcG9sRicvJSdpdGFsaWNHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ2l0YWxpY0YnLUkjbW9HRiQ2LVEjOj1GJy9GOlEnbm9ybWFsRicvJSZmZW5jZUdRJmZhbHNlRicvJSpzZXBhcmF0b3JHRkQvJSlzdHJldGNoeUdGRC8lKnN5bW1ldHJpY0dGRC8lKGxhcmdlb3BHRkQvJS5tb3ZhYmxlbGltaXRzR0ZELyUnYWNjZW50R0ZELyUnbHNwYWNlR1EsMC4yNzc3Nzc4ZW1GJy8lJ3JzcGFjZUdGUy1JJm1mcmFjR0YkNigtRiM2Li1JK211bmRlcm92ZXJHRiQ2Jy1GPTYtUSomUHJvZHVjdDtGJ0ZARkJGRS9GSEY4RkkvRkxGOC9GTkY4Rk8vRlJRJjAuMGVtRicvRlVRLDAuMTY2NjY2N2VtRictRiM2Ji1GLDYlUSJqRidGNkY5LUY9Ni1RIj1GJ0ZARkJGRUZHRklGS0ZNRk9GUUZULUkjbW5HRiQ2JFEiMUYnRkBGQC1GIzYmLUYsNiVRImlGJ0Y2RjktRj02LVEoJm1pbnVzO0YnRkBGQkZFRkdGSUZLRk1GTy9GUlEsMC4yMjIyMjIyZW1GJy9GVUZncEZqb0ZARk8vJSxhY2NlbnR1bmRlckdGREYrLUknbXNwYWNlR0YkNiYvJSdoZWlnaHRHUSYwLjBleEYnLyUmd2lkdGhHUSQ1LjBGJy8lJmRlcHRoR0ZgcS8lKmxpbmVicmVha0dRJWF1dG9GJy1JKG1mZW5jZWRHRiQ2JC1GIzYmLUYsNiVRInRGJ0Y2RjlGY3AtSSVtc3ViR0YkNiUtRiw2JVEieEYnRjZGOS1GIzYkRmRvRkAvJS9zdWJzY3JpcHRzaGlmdEdRIjBGJ0ZARkAtRj02LVEifkYnRkBGQkZFRkdGSUZLRk1GT0Zeby9GVUZfby1GPTYtUScmc2RvdDtGJ0ZARkJGRUZHRklGS0ZNRk9GXm9GX3NGXHMtRmZuNidGaG4tRiM2J0Zkb0Znby1GIzYmRmBwLUY9Ni1RIitGJ0ZARkJGRUZHRklGS0ZNRk9GZnBGaHBGam9GQEYrRkAtRiw2JVEibkYnRjZGOUZPRmlwRitGW3FGaXFGQC1GIzYtRmVuRlxzLUZqcTYkLUYjNiYtRmJyNiVGZHItRiM2JEZgcEZARmlyRmNwRmFyRkBGQEZcc0Zgc0Zcc0Zjc0YrRltxRmF0RkAvJS5saW5ldGhpY2tuZXNzR0ZdcC8lK2Rlbm9tYWxpZ25HUSdjZW50ZXJGJy8lKW51bWFsaWduR0ZddS8lKWJldmVsbGVkR0ZERkBGK0ZARitGQA==pol := product(t-x[j], j = 1 .. i-1)*product(t-x[j], j = i+1 .. n)/(product(x[i]-x[j], j = 1 .. i-1)*product(x[i]-x[j], j = i+1 .. n)):for i from 1 to n do psi[i]:= unapply(pol,t) od; ##unassign('i'): ##c:=y: ##phi:= unapply(sum(y[i]*psi[i](t), i=1..n), t); ##Go to function definition, table of points, choice of approximation: Lagrange, Newton 1, Newton 2, interp, spline, plain, , Taylor,different forms of output, graphical outputs for a given function , table of points.2) Polynomial interpolation - Newton's 1st formulametoda := "NEWTON 1"; k := n; c := array(1 .. n);pol:= proc(i,t) local j, r; r:=c[i];
for j from i-1 by -1 to 1 do r:=c[j]+(t-x[j])*r od; r
end: ##for i from 1 to n do psi[i]:= unapply(pol(i,t), t) od: ##unassign('i'): ##phi:= psi[n]: ##eval(phi);koeficienty();Go to function definition, table of points, choice of approximation: Lagrange, Newton 1, Newton 2, interp, spline, plain, , Taylor,different forms of output, graphical outputs for a given function , table of points.3) Polynomial interpolation - Newton's 2nd formulametoda:="NEWTON 2": k:=n; c:=array(1..n): ##unassign('i','j'): ##for i from 1 to n do psi[i]:= unapply(product(t-x[j], j=1..i-1), t) od; ##unassign('i','j'): ##phi:= unapply(sum(c[i]*psi[i](t), i=1..n), t); ##;koeficienty();4) Polynomial interpolation - standard function interpmetoda:='INTERP': k:=0: ##phi:= unapply(interp(x,y,t),t); ##Go to function definition, table of points, choice of approximation: Lagrange, Newton 1, Newton 2, interp, spline, plain, , Taylor,different forms of output, graphical outputs for a given function , table of points.5) Spline interpolationrad:=3; # order of splinemetoda:=cat("SPLINE (n=", n, ", order=", rad,")"): k:=0: ##phi:= unapply(spline(x,y,t,rad),t); ##Go to function definition, table of points, choice of approximation: Lagrange, Newton 1, Newton 2, interp, spline, plain, , Taylor,different forms of output, graphical outputs for a given function , table of points.6) Plain interpolationChoice of functions psi[i], i=1..nAssigment of powers can be found in Appendix.psi[1]:= 1;psi[2]:= t -> t;psi[3]:= t -> t^2;psi[4]:= t -> sin(t/2)^4;psi[5]:= sin^3;k:=n:Plot of functions psi[i] for plain interpolationplot([psi[i](t) $i=1..k], t=a..b);Computation of interpolating polynomialmetoda:="PLAIN": k:=n; ##c:=array(1..n): ##phi:= unapply(sum(c[i]*psi[i](t), i=1..n), t); ##koeficienty();Go to function definition, table of points, choice of approximation: Lagrange, Newton 1, Newton 2, interp, spline, plain, , Taylor,different forms of output, graphical outputs for a given function , table of points.7) Least Squares Approximation (LSM)Choice of functions psi[i], i=1..k, and their number, k Assigment of powers can be found in Appendix.Plot points. DO NOT SKIP THIS STEP (because of future use)!plot_xy:=plot([[x[i],y[i]] $i=1..n], interval, style=POINT, color=BLUE):
plot_xy;Assignment of functions. k:=3; psi[1]:= 1;psi[2]:= t -> sin((t+1)/2)^4;psi[3]:= sin^3;Plot of functions psi[i] for LSM. DO NOT SKIP THIS STEP (because of future use)!if k>1 then plot([psi[i](t) $i=1..k], t=interval, color=[seq(COLOR(RGB, 1, 1-i/(k-1), 0), i=0..k-1)])
else plot(psi[1](t), t=interval, color=yellow)
end if; plot_psi:=%:Computation of aproximating functionmetoda:="MNC": ##MNCoption:="new"; # choice of library: "new"\342\206\222LinearAlgebra or "old"\342\206\222linlagc:=array(1..k); ##phi:= unapply(sum(c[i]*psi[i](t), i=1..k), t); ##koeficienty();Go to function definition, table of points, choice of approximation: Lagrange, Newton 1, Newton 2, interp, spline, plain, , Taylor,different forms of output, graphical outputs for a given function , table of points.8) Taylor's approximationstred:=3; # center Order:=6; # ordermetoda := cat("TAYLOR (order=",Order, ", center=", stred, ")"); k := 0:rada:=evalf(taylor(f(t), t=stred)); ##phi:=unapply(convert(rada, polynom), t); ##Go to function definition, table of points, choice of approximation: Lagrange, Newton 1, Newton 2, interp, spline, plain, , Taylor,different forms of output, graphical outputs for a given function , table of points.Appendix: Assignment of powers For polynomial aproximation (in plain interpolation or LSM):k:=3:
for i from 1 to k do psi[i]:=unapply(t^(i-1), t) od; i:='i':Go to plain intrepolation, . Different forms of result(Can be customized.)Convert expresion to functional operator:psi:=unapply(expand(phi(t)), t);Resulting interpolation/aproximation:phi(t);Evaluate using floating-point arithmetic:evalf(phi(t));Sort polynomial coefficients:sort(phi(t));Expand and sort:sort(expand(phi(t)));Convert formulae to Horner's form:convert(phi(t), horner, t);Optimization procedure:phi:=codegen[optimize](phi);Testing numerical values:phi(2.5); Go to function definition, table of points, choice of approximation: Lagrange, Newton 1, Newton 2, interp, spline, plain, , Taylor,different forms of output, graphical outputs for a given function , table of points.Graphical representationOption: If you want to change the range of x coordinates, replace a, b by other values.interval:=a..b;a) Given function and interval (you choose points)Graph of interpolation:plot_f:=plot(f(t), t=a..b, color=YELLOW, legend="given function"): plot_phi_l:=plot(phi(t), t=a..b, color=GREEN, title=cat(metoda," n = ", convert(n, string)), legend="interpolation" ): plot_phi:=plot(phi(t), t=a..b, color=GREEN, title=metoda):plot_xy_l:=plot([[x[i],y[i]] $i=1..n], a..b, style=POINT, symbol=POINT, color=BLUE, legend="given points"):plot_xy:=plot([[x[i],y[i]] $i=1..n], a..b, style=POINT, symbol=POINT, color=BLUE):display([plot_f, plot_phi_l, plot_xy_l], title=metoda);Go to function definition, table of points, choice of approximation: Lagrange, Newton 1, Newton 2, interp, spline, plain, , Taylor,different forms of output, graphical outputs for a given function , table of points.Graph of terms of interpolating polynomial (use plot_phi, plot_xy from previous inputs)(uses plot_phi, plot_xy from previous inputs)DO NOT USE for (interp, spline, Taylor)plot_psi:=plot([c[i]*psi[i](t) $i=1..k], t=a..b, -2..3/2, color=[seq(COLOR(RGB, 1, 1-i/(n-1), 0), i=0..n-1)]):display([plot_psi, plot_xy, plot_phi], title=metoda);Alternative recommended for display of consecutive interpolations in NEWTON 1:plot_psi:=plot([psi[i](t) $i=1..k], t=a..b, -2..3/2, color=[seq(COLOR(RGB, 1, 1-i/(n-1), 0), i=0..n-1)]):display([plot_psi, plot_xy, plot_phi], title=metoda);Graph of absolute error(comment out unused commands)without error band epsplot((phi-f)(t), t=a..b, title=cat(metoda," - absolute error"));with error band epseps:=1E-6:plot([(phi-f)(t), -eps, eps], t=interval, -eps..eps, title=cat(metoda," - absolute error"), color=[red, green, green]);Graph of relative error (possibly assign interval for y axis, e.g. -0.1 .. 0.2, ) without error band epsplot((phi/f-1)(t), t=interval, -0.1 .. 0.2, title=cat(metoda," - relative error"));with error band epseps:=1E-6:plot([(phi/f-1)(t), -eps, eps], t=interval, -eps..eps, title=cat(metoda," - relative error"), color=[red, green, green]);Go to function definition, table of points, choice of approximation: Lagrange, Newton 1, Newton 2, interp, spline, plain, , Taylor,different forms of output, graphical outputs for a given function , table of points.JSFHError bound of POLYNOMIAL INTERPOLATION/TAYLOR's APROXIMATIONnth derivative:derf := (D@@n)(f);plot(derf, a..b, title=cat(n, "th derivative"));Bound of nth derivative:Mn := max(evalf(maximize(derf(t), t=a..b)),
-evalf(minimize(derf(t), t=a..b))); ##A) Error bound of POLYNOMIAL INTERPOLATION odhad:= t -> Mn*abs(product(t-x[i], i = 1 .. n))/factorial(n);B) Error bound of TAYLOR's APROXIMATIONodhad:= t -> Mn*abs(t-stred)^n/factorial(n);A and B (together) Error bound and actual error:plot_odhad:=plot(odhad, a..b, title=cat(metoda," - error bound"), color=RED, legend="border on error"):plot_chyba:=plot(abs@(phi-f), a..b, title=cat(metoda," - error bound"), color=GREEN, legend="actual error"):display({plot_odhad, plot_chyba});Go to function definition, table of points, choice of approximation: Lagrange, Newton 1, Newton 2, interp, spline, plain, , Taylor,different forms of output, graphical outputs for a given function , table of points.JSFHGraph of absolute error of differentiationplot(D(phi)-D(f), a..b, title=cat(metoda," - differentiation error"));Go to function definition, table of points, choice of approximation: Lagrange, Newton 1, Newton 2, interp, spline, plain, , Taylor,different forms of output, graphical outputs for a given function , table of points.JSFHGraph of absolute error of integrationplot(int(phi(u), u=x[1]..t)-int(f(u), u=x[1]..t), t=a..b,
title=cat(metoda," - integration error"));Go to function definition, table of points, choice of approximation: Lagrange, Newton 1, Newton 2, interp, spline, plain, , Taylor,different forms of output, graphical outputs for a given function , table of points.JSFHb) Given table of points onlyGraph of aproximation:plot_phi:=plot(phi(t), t=a..b, color=GREEN, title=metoda): plot_xy:=plot([[x[i],y[i]] $i=1..n], a..b, style=POINT, symbol=POINT, color=BLUE):display([plot_phi, plot_xy]);Go to function definition, table of points, choice of approximation: Lagrange, Newton 1, Newton 2, interp, spline, plain, , Taylor,different forms of output, graphical outputs for a given function , table of points.Graf of components of aproximation(uses plot_psi, plot_phi, plot_xy from previous inputs)display([plot_phi, plot_psi, plot_xy]);Graphical representation of LSMinfolevel[Student[LinearAlgebra]] := 1:L:=[seq( [x[i], y[i]], i = 1..n )]:LeastSquaresPlot(L, [t,z], curve=sum(cc[i]*psi[i](t), i=1..k), pointoptions=[symbol=solidcircle, symbolsize=15] );Go to function definition, table of points, choice of approximation: Lagrange, Newton 1, Newton 2, interp, spline, plain, , Taylor,different forms of output, graphical outputs for a given function , table of points.MapletsFor comparison, you may try an interactive Maplet using commands from a new library CurveFitting. CurveFitting[Interactive](x,y);Useful commandsThe following commands can be useful in the solution. An easy example follows (see Help for more details).It is recommended to perform symbolic operations on expressions rather than functions.Expressions can be converted to functions by unapply (see below).Auxiliary expressionsvyraz:=sin(t^3 - 5);uraz:= t^3/10;subs - substitute subexpressions in an expressionsubs(t^3=uraz, vyraz);It can also evaluate an expression at a point. Notice that substitution in function sin does not perform full evaluation and t remains unassigned.subs(t=2., vyraz);
subs(t=2., uraz);
t;eval \342\200\223 evaluation of an expression at a pointeval(vyraz, t=2.);eval(uraz, t=2.);t;Conversion of an expression to a function: unapplyfcevyraz:=unapply(vyraz, t);For a given point, it computes the function value (notice the effect of a decimal point in the argument):fcevyraz(2), fcevyraz(2.);The same construction with an expression creates a meaningless output:vyraz(2), vyraz(2.);piecewise - piecewise-defined functionspom:=piecewise(t>0 and t<2, vyraz, t>=2 and t<5, uraz, 12*t);A function can be defined using unapply ...fce_pom:=unapply(pom, t);
... and we test its correctnessfce_pom(-2), fce_pom(1), fce_pom(1.), fce_pom(3);
In contrast to an expressionpom(2);