{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "````{panels}\n", "Voraussetzungen\n", "^^^\n", "mathematische Grundlagen von\n", "- Newton-Verfahren\n", "- Jacobi-Matrix\n", "- finite Differenzen\n", "---\n", "\n", "Lerninhalte\n", "^^^\n", "- Implementierung und Anwendung des mehrdimensionalen Newton-Verfahrens als selbst erstellte Routine\n", "````" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(newton_aufgaben)=\n", "# Newton-Verfahren\n", "\n", "## Aufgabe 1\n", "\n", "```{image} https://upload.wikimedia.org/wikipedia/commons/thumb/a/a5/A_glass_of_red_wine.jpg/640px-A_glass_of_red_wine.jpg\n", ":width: 500px\n", ":align: center\n", "```\n", "
\"A glass of wine\", Davide Restivo, [CC BY-SA 2.0] via Flickr

\n", "\n", "Gegeben sei ein Weinglas mit dem Querschnitt $y= tan(x^2)$. An welcher Stelle auf der y-Achse muss der Eichstrich für einen Inhalt von $V=0.5$ gezogen werden?\n", "\n", "```{image} images/weinglas_funktion.jpg\n", ":width: 500px\n", ":align: center\n", "```\n", "\n", "```{admonition} Hinweis\n", "1. Um das Volumen auszurechnen, können Sie das Glas auch als Rotationskörper um die x-Achse ansehen. \n", " \n", " $$V=\\pi \\cdot \\int _{a}^{b}(f(x))^{2}\\mathrm {d} x$$\n", "2. Die Stammfunktion von $\\arctan(x)$ ist:\n", "\n", " $$x\\arctan \\left( x \\right) -1/2\\,\\ln \\left( 1+{x}^{2} \\right)$$\n", "```\n", "\n", "Verwenden Sie zur Lösung die Matlab-Funktion `fzero`.\n", "\n", "Aus der Hilfe zu fzero:\n", " \n", "```\n", "FZERO Scalar nonlinear zero finding. \n", "X = FZERO(FUN,X0) tries to find a zero of the function FUN near X0, \n", "if X0 is a scalar. \n", "Examples\n", " FUN can be specified using @:\n", " X = fzero(@sin,3)\n", " returns pi.\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "% SPACE FOR SOLUTION\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Aufgabe 2\n", "\n", "Zum Finden einer Nullstelle implementieren Sie nun das Newton-Verfahren. Angefangen mit einem Startwert $\\mathbf{x}_{0}$ muss in jeder Iteration folgende Verfahrensvorschrift erfüllt werden:\n", "\n", "$$\n", "\\begin{align}\n", " \\mathbf{x}_{i+1} = \\mathbf{x}_i - \\frac{f(\\mathbf{x}_i)}{f^{\\prime}(\\mathbf{x}_i)}, \\quad i=0,1,2,... \\notag\n", "\\end{align}\n", "$$\n", "\n", "a) Programmieren Sie das Newton-Verfahren, so dass es wie `fzero` angewendet werden kann. Da die Ableitung der Funktion noch sehr einfach von Hand zu berechnen ist, kann sie dem Programm übergeben werden.\n", "\n", "```{admonition} Hinweis\n", "Machen Sie sich mit *[Function Handles](https://de.mathworks.com/help/matlab/matlab_prog/creating-a-function-handle.html)* und *[Anonymous Functions](https://www.mathworks.com/help/matlab/matlab_prog/anonymous-functions.html)* in Matlab vertraut.\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%file simple_newton.m\n", "function x = simple_newton(func,deriv,x0,tol,maxit)\n", "% z = simple_newton(f,df,x0,tol,maxit) solves the nonlinear system 0=func(x)\n", "%\n", "% inputs:\n", "% func a handle to the nonlinear function\n", "% deriv a handle to the derivative of the nonlinear function\n", "% x0 initial guess for the Newton method\n", "% atol absolute tolerance\n", "% maxit maximum number of Newton iterations\n", "\n", "% YOUR CODE HERE" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Das erstellte Programm lässt sich folgendermaßen testen:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "moxunit_runtests test_simple_newton.m" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "b) Lösen Sie Aufgabe 1 noch einmal mit Hilfe Ihres soeben programmierten Newton-Verfahrens und schauen Sie, ob die Ergebnisse übereinstimmen." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "% SPACE FOR SOLUTION\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Aufgabe 3\n", "\n", "Gegeben sei das nichtlineare System\n", "\n", "$$\n", "\\begin{align}\n", "5x_1^2-x_2^2 &= 0 \\\\\n", "x_1 x_2-\\frac{1}{4}(sin(x_1)+cos(x_2)) &= -\\frac{1}{4}\n", "\\end{align}\n", "$$\n", "\n", "und eine Lösung für $x$ liegt in der Nähe von \n", "\n", "$$\n", "\\begin{align}\n", "\\left( \\begin{array}{c} x_1\\\\x_2 \\end{array} \\right) &=\n", "\\left( \\begin{array}{c} 1/8\\\\1/4 \\end{array} \\right)\n", "\\end{align}.\n", "$$\n", "\n", "a) Bestimmen Sie die Jacobimatrix des Systems.\n", "\n", "Um dem Newton-Verfahren nicht jedes mal eine Jacobi-Matrix übergeben zu müssen, kann diese auch mittels finiter Differenzen approximiert werden. Das hat den Vorteil, dass so auch komplexe Systeme einfach gelöst werden können. Aufgrund der Approximierung verlieren wir aber an Genauigkeit. \n", "\n", "b) Schreiben Sie eine Matlab-Funktion, die die Jacobimatrix $J_F(\\mathbf{x})$ für eine beliebige Funktion $F:\\mathbb{R}^n \\to \\mathbb{R}^m$ mit Differenzenquotienten (*\"finiten Differenzen\"*) approximiert. Dabei soll die Eingabe `F` ein *[Function Handle](https://de.mathworks.com/help/matlab/matlab_prog/creating-a-function-handle.html)* einer beliebigen Funktion sein, beispielsweise `@sin` oder in diesem Fall `@(x) [5*x(1)^2-x(2)^2; x(1)*x(2)-1/4*(sin(x(1))+cos(x(2)))+1/4]`. `x` ist hierbei der Vektor, an dem die Matrix ausgewertet werden soll und wird als `[1/8; 1/4]` übergeben." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```{admonition} Hinweis\n", "- Sie können Ihr Programm in der Aufgabe {ref}`bierschaum_newton` wiederverwenden. Speichern Sie Ihre Lösung also bestenfalls nach der Fertigstellung ab.\n", "- Implementieren Sie eine Schleife über die Spalten der Jacobi-Matrix, pro Spalte müssen Sie einen Differenzenquotienten bilden.\n", "- Ihr Programm sollte mit $n+1$ Funktionsaufrufen von `F` auskommen.\n", "- Achten Sie darauf, dass alle Funktionen $F:\\mathbb{R}^n \\to \\mathbb{R}^m$ auf die Sie die `jacobian` Funktion anwenden *Spalten*vektoren auf *Spalten*vektoren abbilden.\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%file jacobian.m\n", "function J = jacobian(F,x)\n", "% J = jacobian(F,x) returns the (m x n) Jacobian matrix of F evaluated at x\n", "%\n", "% | dF1/dx1 ... dF1/dxn |\n", "% | . . |\n", "% | dFm/dx1 ... dFm/dxn |\n", "%\n", "% It uses finite difference approximations. x must be a (n,1)-column vector and F must be a function\n", "% taking an (n,1)-vector as an input. m is deduced from F.\n", "\n", "% YOUR CODE HERE" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Prüfen Sie nun mit Hilfe Ihrer in a) errechneten Jacobimatrix die Ergebnisse für verschiedene x. Sie können Ihr Programm des Weiteren auch mittels eines Unit Tests überprüfen:\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "moxunit_runtests test_jacobian.m" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "b) Erweitern Sie mit Hilfe des zuvor erstellten Programms zur Auswertung der Jacobimatrix Ihr Programm für das Newton-Verfahren, sodass es nun das oben genannte Gleichungssystem lösen kann." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```{admonition} Hinweis\n", "Sie können Ihr Programm in der Aufgabe {ref}`bierschaum_newton` wiederverwenden. Speichern Sie Ihre Lösung also bestenfalls nach der Fertigstellung ab.\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "%%file newton.m\n", "function z = newton(func,z0,tol,maxit)\n", "% z = newton(F,z0,tol,maxit) solves the nonlinear system 0=func(z)\n", "%\n", "% inputs:\n", "% func a handle to the nonlinear function\n", "% z0 initial guess for the Newton method\n", "% atol absolute tolerance\n", "% maxit maximum number of Newton iterations\n", "\n", "% YOUR CODE HERE" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Auch dieses Programm lässt sich mittels eines Unit Tests überprüfen:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "moxunit_runtests test_newton.m" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Ermitteln Sie nun die Lösung des oben genannten Gleichungssystems:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "% solve nonlinear function for x\n", "x = newton(@(x) [5*x(1)^2-x(2)^2; x(1)*x(2)-1/4*(sin(x(1))+cos(x(2)))+1/4],[1/8; 1/4],eps,100)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Zur Kontrolle berechnen wir das Gleichungssystem nochmal rückwärts. \n", "\n", "Erinnerung: \n", "\n", "$$\n", "\\begin{align}\n", "\\left( \\begin{array}{c} y_1\\\\y_2 \\end{array} \\right) &=\n", "\\left( \\begin{array}{c} 0\\\\-1/4 \\end{array} \\right)\n", "\\end{align}\n", "$$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "% check if solution solves the nonlinear system\n", "y = [5*x(1)^2-x(2)^2; x(1)*x(2)-1/4*(sin(x(1))+cos(x(2)))]\n", "if max(abs(y-[0;-1/4])) > eps\n", " disp(['Incorrect solution with a tolerance of |', num2str(max(abs(y-[0;-1/4]))), '| > eps']) \n", "else\n", " disp(['Correct solution with a tolerance of |', num2str(max(abs(y-[0;-1/4]))), '| <= eps'])\n", "end" ] } ], "metadata": { "kernelspec": { "display_name": "Octave", "language": "octave", "name": "octave" }, "language_info": { "file_extension": ".m", "help_links": [ { "text": "GNU Octave", "url": "https://www.gnu.org/software/octave/support.html" }, { "text": "Octave Kernel", "url": "https://github.com/Calysto/octave_kernel" }, { "text": "MetaKernel Magics", "url": "https://metakernel.readthedocs.io/en/latest/source/README.html" } ], "mimetype": "text/x-octave", "name": "octave", "version": "5.2.0" } }, "nbformat": 4, "nbformat_minor": 4 }