Skip to content
Permalink
f88580cf1a
Switch branches/tags

Name already in use

A tag already exists with the provided branch name. Many Git commands accept both tag and branch names, so creating this branch may cause unexpected behavior. Are you sure you want to create this branch?

ME3255F2017/17_integrals_and_derivatives/17_integrals.ipynb

Go to file
 
 
Cannot retrieve contributors at this time
543 lines (543 sloc) 21.9 KB
Raw Blame
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"attributes": {
"classes": [],
"id": "",
"n": "1"
}
},
"outputs": [],
"source": [
"setdefaults"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"attributes": {
"classes": [],
"id": "",
"n": "2"
}
},
"outputs": [],
"source": [
"%plot --format svg"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Numerical Integrals\n",
"\n",
"# Integrals in practice\n",
"\n",
"### Example: Compare toughness of Stainless steel to Structural steel\n",
"\n",
"![Stress-strain plot of steel](steel_psi.jpg)\n",
"\n",
"### Step 1 - G3Data to get points\n",
"\n",
"Use the plot shown to determine the toughness of stainless steel and the\n",
"toughness of structural steel."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"attributes": {
"classes": [],
"id": "",
"n": "9"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"toughness of structural steel is 10.2 psi\n",
"toughness of stainless steel is 18.6 psi\n"
]
}
],
"source": [
"fe_c=load('structural_steel_psi.jpg.dat');\n",
"fe_cr =load('stainless_steel_psi.jpg.dat');\n",
"\n",
"fe_c_toughness=trapz(fe_c(:,1),fe_c(:,2));\n",
"fe_cr_toughness=trapz(fe_cr(:,1),fe_cr(:,2));\n",
"\n",
"fprintf('toughness of structural steel is %1.1f psi\\n',fe_c_toughness)\n",
"fprintf('toughness of stainless steel is %1.1f psi',fe_cr_toughness)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Gauss Quadrature (for functions)\n",
"\n",
"Evaluating an integral, we assumed a polynomial form for each Newton-Cotes\n",
"approximation.\n",
"\n",
"If we can evaluate the function at any point, it makes more sense to choose\n",
"points more wisely rather than just using endpoints\n",
"\n",
"![trapezoidal example](trap_example.png)\n",
"\n",
"Set up two unknown constants, $c_{0}$ and $x_{0}$ and determine a *wise* place\n",
"to evaluate f(x) such that\n",
"\n",
"$I=c_{0}f(x_{0})$\n",
"\n",
"and I is exact for polynomial of n=0, 1\n",
"\n",
"$\\int_{a}^{b}1dx=b-a=c_{0}$\n",
"\n",
"$\\int_{a}^{b}xdx=\\frac{b^2-a^2}{2}=c_{0}x_{0}$\n",
"\n",
"so $c_{0}=b-a$ and $x_{0}=\\frac{b+a}{2}$\n",
"\n",
"$I=\\int_{a}^{b}f(x)dx \\approx (b-a)f\\left(\\frac{b+a}{2}\\right)$"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"attributes": {
"classes": [],
"id": "",
"n": "10"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"f1 =\n",
"\n",
"@(x) x + 1\n",
"\n",
"f2 =\n",
"\n",
"@(x) 1 / 2 * x .^ 2 + x + 1\n",
"\n",
"f3 =\n",
"\n",
"@(x) 1 / 6 * x .^ 3 + 1 / 2 * x .^ 2 + x\n",
"\n"
]
},
{
"data": {
"image/svg+xml": [
"<svg height=\"420px\" viewBox=\"0 0 560 420\" width=\"560px\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
"\n",
"<title>Gnuplot</title>\n",
"<desc>Produced by GNUPLOT 5.0 patchlevel 3 </desc>\n",
"\n",
"<g id=\"gnuplot_canvas\">\n",
"\n",
"<rect fill=\"none\" height=\"420\" width=\"560\" x=\"0\" y=\"0\"/>\n",
"<defs>\n",
"\n",
"\t<circle id=\"gpDot\" r=\"0.5\" stroke-width=\"0.5\"/>\n",
"\t<path d=\"M-1,0 h2 M0,-1 v2\" id=\"gpPt0\" stroke=\"currentColor\" stroke-width=\"0.222\"/>\n",
"\t<path d=\"M-1,-1 L1,1 M1,-1 L-1,1\" id=\"gpPt1\" stroke=\"currentColor\" stroke-width=\"0.222\"/>\n",
"\t<path d=\"M-1,0 L1,0 M0,-1 L0,1 M-1,-1 L1,1 M-1,1 L1,-1\" id=\"gpPt2\" stroke=\"currentColor\" stroke-width=\"0.222\"/>\n",
"\t<rect height=\"2\" id=\"gpPt3\" stroke=\"currentColor\" stroke-width=\"0.222\" width=\"2\" x=\"-1\" y=\"-1\"/>\n",
"\t<rect fill=\"currentColor\" height=\"2\" id=\"gpPt4\" stroke=\"currentColor\" stroke-width=\"0.222\" width=\"2\" x=\"-1\" y=\"-1\"/>\n",
"\t<circle cx=\"0\" cy=\"0\" id=\"gpPt5\" r=\"1\" stroke=\"currentColor\" stroke-width=\"0.222\"/>\n",
"\t<use fill=\"currentColor\" id=\"gpPt6\" stroke=\"none\" xlink:href=\"#gpPt5\"/>\n",
"\t<path d=\"M0,-1.33 L-1.33,0.67 L1.33,0.67 z\" id=\"gpPt7\" stroke=\"currentColor\" stroke-width=\"0.222\"/>\n",
"\t<use fill=\"currentColor\" id=\"gpPt8\" stroke=\"none\" xlink:href=\"#gpPt7\"/>\n",
"\t<use id=\"gpPt9\" stroke=\"currentColor\" transform=\"rotate(180)\" xlink:href=\"#gpPt7\"/>\n",
"\t<use fill=\"currentColor\" id=\"gpPt10\" stroke=\"none\" xlink:href=\"#gpPt9\"/>\n",
"\t<use id=\"gpPt11\" stroke=\"currentColor\" transform=\"rotate(45)\" xlink:href=\"#gpPt3\"/>\n",
"\t<use fill=\"currentColor\" id=\"gpPt12\" stroke=\"none\" xlink:href=\"#gpPt11\"/>\n",
"\t<path d=\"M0,1.330 L1.265,0.411 L0.782,-1.067 L-0.782,-1.076 L-1.265,0.411 z\" id=\"gpPt13\" stroke=\"currentColor\" stroke-width=\"0.222\"/>\n",
"\t<use fill=\"currentColor\" id=\"gpPt14\" stroke=\"none\" xlink:href=\"#gpPt13\"/>\n",
"\t<filter filterUnits=\"objectBoundingBox\" height=\"1\" id=\"textbox\" width=\"1\" x=\"0\" y=\"0\">\n",
"\t <feFlood flood-color=\"white\" flood-opacity=\"1\" result=\"bgnd\"/>\n",
"\t <feComposite in=\"SourceGraphic\" in2=\"bgnd\" operator=\"atop\"/>\n",
"\t</filter>\n",
"\t<filter filterUnits=\"objectBoundingBox\" height=\"1\" id=\"greybox\" width=\"1\" x=\"0\" y=\"0\">\n",
"\t <feFlood flood-color=\"lightgrey\" flood-opacity=\"1\" result=\"grey\"/>\n",
"\t <feComposite in=\"SourceGraphic\" in2=\"grey\" operator=\"atop\"/>\n",
"\t</filter>\n",
"</defs>\n",
"<g color=\"white\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"1.00\">\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"1.00\">\n",
"\t<g shape-rendering=\"crispEdges\" stroke=\"none\">\n",
"\t\t<polygon fill=\"rgb(255, 255, 255)\" points=\"62.2,384.0 534.9,384.0 534.9,16.8 62.2,16.8 \"/>\n",
"\t</g>\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"rgb(255, 255, 255)\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"\t<path d=\"M62.2,384.0 L74.7,384.0 M535.0,384.0 L522.5,384.0 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"end\" transform=\"translate(53.9,390.0)\">\n",
"\t\t<text><tspan font-family=\"{}\">-1000</tspan></text>\n",
"\t</g>\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"\t<path d=\"M62.2,310.5 L74.7,310.5 M535.0,310.5 L522.5,310.5 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"end\" transform=\"translate(53.9,316.5)\">\n",
"\t\t<text><tspan font-family=\"{}\">-500</tspan></text>\n",
"\t</g>\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"\t<path d=\"M62.2,237.1 L74.7,237.1 M535.0,237.1 L522.5,237.1 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"end\" transform=\"translate(53.9,243.1)\">\n",
"\t\t<text><tspan font-family=\"{}\">0</tspan></text>\n",
"\t</g>\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"\t<path d=\"M62.2,163.6 L74.7,163.6 M535.0,163.6 L522.5,163.6 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"end\" transform=\"translate(53.9,169.6)\">\n",
"\t\t<text><tspan font-family=\"{}\">500</tspan></text>\n",
"\t</g>\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"\t<path d=\"M62.2,90.2 L74.7,90.2 M535.0,90.2 L522.5,90.2 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"end\" transform=\"translate(53.9,96.2)\">\n",
"\t\t<text><tspan font-family=\"{}\">1000</tspan></text>\n",
"\t</g>\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"\t<path d=\"M62.2,16.7 L74.7,16.7 M535.0,16.7 L522.5,16.7 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"end\" transform=\"translate(53.9,22.7)\">\n",
"\t\t<text><tspan font-family=\"{}\">1500</tspan></text>\n",
"\t</g>\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"\t<path d=\"M62.2,384.0 L62.2,371.5 M62.2,16.7 L62.2,29.2 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"middle\" transform=\"translate(62.2,408.0)\">\n",
"\t\t<text><tspan font-family=\"{}\">-20</tspan></text>\n",
"\t</g>\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"\t<path d=\"M121.3,384.0 L121.3,371.5 M121.3,16.7 L121.3,29.2 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"middle\" transform=\"translate(121.3,408.0)\">\n",
"\t\t<text><tspan font-family=\"{}\">-15</tspan></text>\n",
"\t</g>\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"\t<path d=\"M180.4,384.0 L180.4,371.5 M180.4,16.7 L180.4,29.2 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"middle\" transform=\"translate(180.4,408.0)\">\n",
"\t\t<text><tspan font-family=\"{}\">-10</tspan></text>\n",
"\t</g>\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"\t<path d=\"M239.5,384.0 L239.5,371.5 M239.5,16.7 L239.5,29.2 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"middle\" transform=\"translate(239.5,408.0)\">\n",
"\t\t<text><tspan font-family=\"{}\">-5</tspan></text>\n",
"\t</g>\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"\t<path d=\"M298.6,384.0 L298.6,371.5 M298.6,16.7 L298.6,29.2 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"middle\" transform=\"translate(298.6,408.0)\">\n",
"\t\t<text><tspan font-family=\"{}\">0</tspan></text>\n",
"\t</g>\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"\t<path d=\"M357.7,384.0 L357.7,371.5 M357.7,16.7 L357.7,29.2 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"middle\" transform=\"translate(357.7,408.0)\">\n",
"\t\t<text><tspan font-family=\"{}\">5</tspan></text>\n",
"\t</g>\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"\t<path d=\"M416.8,384.0 L416.8,371.5 M416.8,16.7 L416.8,29.2 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"middle\" transform=\"translate(416.8,408.0)\">\n",
"\t\t<text><tspan font-family=\"{}\">10</tspan></text>\n",
"\t</g>\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"\t<path d=\"M475.9,384.0 L475.9,371.5 M475.9,16.7 L475.9,29.2 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"middle\" transform=\"translate(475.9,408.0)\">\n",
"\t\t<text><tspan font-family=\"{}\">15</tspan></text>\n",
"\t</g>\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"\t<path d=\"M535.0,384.0 L535.0,371.5 M535.0,16.7 L535.0,29.2 \" stroke=\"black\"/>\t<g fill=\"rgb(0,0,0)\" font-family=\"{}\" font-size=\"16.00\" stroke=\"none\" text-anchor=\"middle\" transform=\"translate(535.0,408.0)\">\n",
"\t\t<text><tspan font-family=\"{}\">20</tspan></text>\n",
"\t</g>\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"\t<path d=\"M62.2,16.7 L62.2,384.0 L535.0,384.0 L535.0,16.7 L62.2,16.7 Z \" stroke=\"black\"/></g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"</g>\n",
"\t<g id=\"gnuplot_plot_1a\"><title>gnuplot_plot_1a</title>\n",
"<g color=\"white\" fill=\"none\" stroke=\"black\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"3.00\">\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"3.00\">\n",
"\t<path d=\"M85.8,358.7 L90.1,351.1 L94.4,343.9 L98.7,336.9 L103.0,330.3 L107.3,324.0 L111.6,317.9 L115.9,312.2 L120.2,306.7 L124.5,301.5 L128.8,296.6 L133.1,291.9 L137.4,287.5 L141.7,283.3 L146.0,279.4 L150.3,275.7 L154.6,272.2 L158.9,269.0 L163.2,265.9 L167.5,263.1 L171.8,260.4 L176.1,258.0 L180.4,255.7 L184.7,253.6 L189.0,251.6 L193.3,249.9 L197.6,248.3 L201.9,246.8 L206.2,245.4 L210.5,244.2 L214.8,243.2 L219.1,242.2 L223.4,241.4 L227.7,240.6 L232.0,240.0 L236.3,239.4 L240.6,238.9 L244.9,238.5 L249.2,238.2 L253.5,237.9 L257.8,237.7 L262.1,237.6 L266.4,237.4 L270.7,237.3 L275.0,237.3 L279.3,237.2 L283.6,237.2 L287.9,237.2 L292.2,237.1 L296.5,237.1 L300.7,237.1 L305.0,237.0 L309.3,236.9 L313.6,236.7 L317.9,236.5 L322.2,236.3 L326.5,236.0 L330.8,235.6 L335.1,235.2 L339.4,234.7 L343.7,234.1 L348.0,233.4 L352.3,232.6 L356.6,231.7 L360.9,230.7 L365.2,229.5 L369.5,228.3 L373.8,226.9 L378.1,225.3 L382.4,223.6 L386.7,221.8 L391.0,219.7 L395.3,217.5 L399.6,215.2 L403.9,212.6 L408.2,209.9 L412.5,206.9 L416.8,203.8 L421.1,200.4 L425.4,196.8 L429.7,193.0 L434.0,189.0 L438.3,184.7 L442.6,180.1 L446.9,175.3 L451.2,170.3 L455.5,164.9 L459.8,159.3 L464.1,153.4 L468.4,147.2 L472.7,140.8 L477.0,134.0 L481.3,126.9 L485.6,119.5 L489.9,111.7 L494.2,103.6 L498.5,95.2 L502.8,86.4 L507.1,77.3 L511.4,67.8 \" stroke=\"rgb( 0, 0, 255)\"/></g>\n",
"\t</g>\n",
"<g color=\"white\" fill=\"none\" stroke=\"rgb( 0, 0, 255)\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"2.00\">\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"2.00\">\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"black\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"</g>\n",
"<g color=\"black\" fill=\"none\" stroke=\"currentColor\" stroke-linecap=\"butt\" stroke-linejoin=\"miter\" stroke-width=\"0.50\">\n",
"</g>\n",
"</g>\n",
"</svg>"
],
"text/plain": [
"<IPython.core.display.SVG object>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"f1=@(x) x+1\n",
"f2=@(x) 1/2*x.^2+x+1\n",
"f3=@(x) 1/6*x.^3+1/2*x.^2+x\n",
"plot(linspace(-18,18),f3(linspace(-18,18)))"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"attributes": {
"classes": [],
"id": "",
"n": "11"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"integral of f1 from 2 to 3 = 3.500000\n",
"integral of f1 from 2 to 3 ~ 3.500000\n",
"integral of f2 from 2 to 3 = 6.666667\n",
"integral of f2 from 2 to 3 ~ 6.625000\n"
]
}
],
"source": [
"fprintf('integral of f1 from 2 to 3 = %f',f2(3)-f2(2))\n",
"fprintf('integral of f1 from 2 to 3 ~ %f',(3-2)*f1(3/2+2/2))\n",
"\n",
"fprintf('integral of f2 from 2 to 3 = %f',f3(3)-f3(2))\n",
"fprintf('integral of f2 from 2 to 3 ~ %f',(3-2)*f2(3/2+2/2))\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This process is called **Gauss Quadrature**. Usually, the bounds are fixed at -1\n",
"and 1 instead of a and b\n",
"\n",
"$I=c_{0}f(x_{0})$\n",
"\n",
"and I is exact for polynomial of n=0, 1\n",
"\n",
"$\\int_{-1}^{1}1dx=b-a=c_{0}$\n",
"\n",
"$\\int_{-1}^{1}xdx=\\frac{1^2-(-1)^2}{2}=c_{0}x_{0}$\n",
"\n",
"so $c_{0}=2$ and $x_{0}=0$\n",
"\n",
"$I=\\int_{-1}^{1}f(x)dx \\approx 2f\\left(0\\right)$\n",
"\n",
"Now, integrals can be performed with a change of variable\n",
"\n",
"a=2\n",
"\n",
"b=3\n",
"\n",
"x= 2 to 3\n",
"\n",
"or $x_{d}=$ -1 to 1\n",
"\n",
"$x=a_{1}+a_{2}x_{d}$\n",
"\n",
"at $x_{d}=-1$, x=a\n",
"\n",
"at $x_{d}=1$, x=b\n",
"\n",
"so\n",
"\n",
"$x=\\frac{(b+a) +(b-a)x_{d}}{2}$\n",
"\n",
"$dx=\\frac{b-a}{2}dx_{d}$\n",
"\n",
"$\\int_{2}^{3}x+1dx=\\int_{-1}^{1}\\left(\\frac{(2+3) +(3-2)x_{d}}{2}\n",
"+1\\right)\n",
"\\frac{3-2}{2}dx_{d}$\n",
"\n",
"$\\int_{2}^{3}x+1dx=\\int_{-1}^{1}\\left(\\frac{5 +x_{d}}{2}\n",
"+1\\right)\n",
"\\frac{3-2}{2}dx_{d}$\n",
"\n",
"$\\int_{2}^{3}x+1dx=\\int_{-1}^{1}\\left(\\frac{7}{4}+\\frac{1}{4}x_{d}\\right)dx_{d}=\n",
"2\\frac{7}{4}=3.5$"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"attributes": {
"classes": [],
"id": "",
"n": "12"
}
},
"outputs": [],
"source": [
"function I=gauss_1pt(func,a,b)\n",
" % Gauss quadrature using single point\n",
" % exact for n<1 polynomials\n",
" c0=2;\n",
" xd=0;\n",
" dx=(b-a)/2;\n",
" x=(b+a)/2+(b-a)/2*xd;\n",
" I=func(x).*dx*c0;\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"attributes": {
"classes": [],
"id": "",
"n": "14"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"ans = 3.5000\r\n"
]
}
],
"source": [
"gauss_1pt(f1,2,3)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## General Gauss weights and points\n",
"\n",
"![Gauss quadrature table](gauss_weights.png)\n",
"\n",
"### If you need to evaluate an integral, to increase accuracy, increase number\n",
"of Gauss points\n",
"\n",
"### Adaptive Quadrature\n",
"\n",
"Matlab/Octave built-in functions use two types of adaptive quadrature to\n",
"increase accuracy of integrals of functions.\n",
"\n",
"1. `quad`: Simpson quadrature good for nonsmooth functions\n",
"\n",
"2. `quadl`: Lobatto quadrature good for smooth functions"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"attributes": {
"classes": [
"matlab"
],
"id": ""
}
},
"outputs": [],
"source": [
"q = quad(fun, a, b, tol, trace, p1, p2, …)\n",
"fun : function to be integrates\n",
"a, b: integration bounds\n",
"tol: desired absolute tolerance (default: 10-6)\n",
"trace: flag to display details or not\n",
"p1, p2, …: extra parameters for fun\n",
"quadl has the same arguments"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"attributes": {
"classes": [],
"id": "",
"n": "15"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"ans = 6.6667\n",
"ans = 6.6667\n",
"ans = 6.6667\n"
]
}
],
"source": [
"% integral of quadratic\n",
"quad(f2,2,3)\n",
"quadl(f2,2,3)\n",
"f3(3)-f3(2)"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"attributes": {
"classes": [],
"id": "",
"n": "16"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"f_c =\n",
"\n",
"@(x) cosh (x / 30) + exp (-10 * x) .* erf (x)\n",
"\n",
"ans = 2.0048\n",
"ans = 2.0048\n"
]
}
],
"source": [
"f_c=@(x) cosh(x/30)+exp(-10*x).*erf(x)\n",
"\n",
"quad(f_c,1,3)\n",
"quadl(f_c,1,3)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Thanks"
]
}
],
"metadata": {},
"nbformat": 4,
"nbformat_minor": 2
}