{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Lecture 3 Part 2 - Non Linear Equations - In-Class Exercises\n", "\n", "\n", "## Exercise 1 - Relaxation method\n", "\n", "Solve the equations a) $x=2-e^{-x}$ and b) $x=e^{1-x^2}$. \n", "\n", "1. Write a short program that iterates equations a) and b).\n", "2. Start iterating a) from $x=1.0$ and b) from $x=0.5$.\n", "3. At every step, take the difference to the value at the previous step and plot this difference as a function of iteration number.\n", "\n", "\n", "## Talking points\n", "\n", "1. What do you observe?\n", "2. How many iterations do you need for 1E-3 accuracy?\n", "3. What is happening in case b)?\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Example Ferromagnetism\n", "\n", "In the mean-field theory of ferromagnetism, the strength $M$ of magnetization of a ferromagnetic material like iron depends on temperature $T$ according to $$M=\\mu \\tanh{\\frac{JM}{k_B T}},$$ where $\\mu$ is the magnetic moment, $J$ is a coupling constant, and $k_B$ is Boltzmann's constant. To simplify things, let us make the substitutions $m=M/\\mu$ and $C=\\mu J/k_B$ so that $$m=\\tanh{\\frac{Cm}{T}}.$$ It is clear that this equation has the trivial solution $m=0$, which implies that a material that is not magnetized at all. But are there solutions with $m\\ne 0$? Here is a program to find the solutions and to make a plot as a function of temperature for $C=1$: " ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEQCAYAAACeDyIUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAsO0lEQVR4nO3deXxU5dn/8c+VBULCFvZFICAIAiKVaMENqrbiBvbXutRWrXWvttrtqT5tqVpb28fWpdXHirVq1T5al4q2VdxREdqC1WoUEFlFQHbCTuD6/XFOdDpOMpnJZE5m5vt+veY1yTnnnnPN8ciV+9ybuTsiIiKpKoo6ABERyU1KICIikhYlEBERSYsSiIiIpEUJRERE0lISdQDZ1K1bN6+qqoo6DBGRnDF37ty17t490b6CSiBVVVXMmTMn6jBERHKGmS1taJ8eYYmISFqUQEREJC1KICIikpasJxAz28fMfmNms8xsm5m5mVU1sWyZmV1vZivNbHv4GUe2cMgiIpJAFDWQwcCpwAbg5RTL3gmcD0wBTgRWAtPNbHQmAxQRkeSi6IX1krv3BDCz84DPNaWQmR0InAF8zd3vCrfNAGqAa4BJLROuiIgkkvUaiLvvTbPoJGA38GDMZ9UBDwDHmlnbDIQnIiJNlEvjQEYAi919W9z2GqANwaOxmpY48a+fe5e6PXvBDAPMgu2GYQbhr8HP9Tvrf090TLjt4+Msbn+47aPz8NG5AYqLjLLSItqVFtOuTUnwXlpMuzZFtG9bSmVFKW1LilvgSoiIfCyXEkgXgnaTeOtj9n+CmV0AXADQv3//tE782xnvsW3XnrTKRqV92xIqK0rpUtGWbhVt2KeyHftUln/0XtWtnA5lpVGHKSI5LJcSiAGJVr+yBNs+4u5TgakA1dXVaa2e9fY1E+M/E/ePg3F3HKhfm8txYtfpCo71mP0fl6nfT3jMx8c3/Nl798KO3XvYXv/aFb5272Hzjt1s2LqL9Vt3s37rTtZv282Kjdv5x+L11O6s+4/v0bdzO4b16sB+vTowvHdHxgyopE/ndulcIhEpQLmUQNYDiaoQlTH7syL28VK4JVunbpZN23fz/oZtLF+/nffWbGHeqloWrKplxoI11O0NMlSfTmWMqerCwVWVHDmkO1XdKiKOWkRaq1xKIDXA582sPK4dZDiwC1gYTVi5o1O7Ujq168SIPp3+Y/uuur3MW7WZuUs3MGfpBv65eD1PvPEBAFVdy5kwtAdHDevBuH27UlqssaciErAo10QPu/HeAQx09yVJjh0N/Av4qrvfE24rAd4EFrr7ScnOV11d7ZpMMTl3Z9n6bbw4fw0vzP+QWe+tY2fdXirLSzlhVG8mHdiX6gGVFBXlRs1LRNJnZnPdvTrRvkhqIGb2xfDHMeH7cWa2Bljj7jPMbADwHnCNu18D4O6vm9mDwE1mVgosBi4GBgJfzu43yG9mxoCuFZx9aAVnH1rFjt17eGnBGh5/4wMenvs+981eRt/O7Tj94H6cdkg/enQoizpkEYlAJDUQM2vopDPcfUI4tcli4Gp3vyqmXDvgpwQDCjsDbwDfd/cXm3Je1UCab+vOOp59ZzUPzXmfVxaupaTIOHZkL845tIrqqoQd4UQkhzVWA4n0EVa2KYFk1qI1W7j/78t4aM5yNu+o45CBXbjkM4M5cki3/xgPIyK5SwkkpATSMrbv2sMD/1zG1JcWsXLTDkb27ch3PjuUCUO7K5GI5DglkJASSMvaVbeXx/61gltfXMjSddsYN6grVx4/jFH7dI46NBFJU2MJRH0yJWPalBRx6sH9eOZb47l60ggWrK5l0i0z+eb//YvVm3dEHZ6IZJgSiGRcm5Iizj60ihe/N4FLPzOY6TWrOPpXM7h75mL27C2cGq9IvlMCkRbToayU7x47lKe/dSSf6t+Zq554m5NvnclbKzZFHZqIZIASiLS4AV0r+MPXDuHXX/oUKzft4ORbZ3LrCwtVGxHJcUogkhVmxqQD+/Dst49k4sheXD99PqfePotl6+Jn5xeRXKEEIlnVubwNt5xxEDefPpoFq2s57uaXeOxfK6IOS0TSoAQikZg8ui/TLz+SEX06cfmDrzNl2lvsqkt3sUoRiYISiESmT+d23H/+p7ngyEH8YdZSTr19Fh9s3B51WCLSREogEqnS4iL++/j9+d8vH8S74biRN5ZvjDosEWkCJRBpFY4/oDePXXIYZaVFnDZ1Fk++uTLqkEQkCSUQaTWG9OzAY5ccxv69O3Lx/a9x+4z3KKSpdkRyjRKItCrd2rfl/84fywmjenPdk/O47sl5SiIirVQuLWkrBaKstJjfnP4pula0YepLi6jdUce1J4+kWCsgirQqSiDSKhUVGVdPGkHHslJueWEhW3bWccOpB2pNdpFWRAlEWi0z47vHDqVDWQnXPTmP3XV7ueWMT1GiJCLSKuj/RGn1Lhy/Lz8+aThP1aziW396Q3NoibQSqoFITjjnsIHsrNvLz5+cR5viIq7/4iiK1CYiEiklEMkZF43fl111e7nhmQWUlRZx7ckjtWSuSISUQCSnfOOowWzbtYffzniPXh3L+MbRQ6IOSaRgKYFITjEzvj9xKB/W7uBXzyygV6cyTqnuF3VYIgVJCURyjpnx8/83ijW1O7ni0Tfp3qEtE4b2iDoskYKjXliSk9qUFHHbV8YwtGcHvn7/a7yzcnPUIYkUHCUQyVnt25Zw1zkH06GshPP/MIf1W3dFHZJIQVECkZzWs2MZt59ZzYe1O7nk/tfYvUeLUolkixKI5LzR/Tpz3ecPYNaidfz0r+9EHY5Iwch6AjGzfmb2sJltMrPNZvaomfVvYtn+ZnaPmS0zs21mtsDMrjWzipaOW1q3L4zZh3MPH8jdry7hkbnvRx2OSEHIai8sMysHngd2AmcDDlwLvGBmo9x9ayNlK4BngVLgR8Ay4GDgamAIcFrLRi+t3ZXHDePtDzbzw8feYtQ+nRjSs0PUIYnktWzXQM4HBgEnu/tj7j4NmAQMAC5MUvYwgkRxobvf4+4vuPv/ADcDXwiTkxSwkuIibj59NOVtirn0j/9i+649UYckkteynUAmAbPdfWH9BndfDMwEJicp2yZ8j++vuZHge2hOC6FHxzJuPG0081fXcvUTNVGHI5LXsp1ARgBvJdheAwxPUvZZ4F3gF2Y23Mzam9lRwGXAbxt7/CWF5cj9uvP1CfvywD+XM+31FVGHI5K3sp1AugAbEmxfD1Q2VtDddwCHE8RcA9QCzwF/AS5tqJyZXWBmc8xszpo1a9KNW3LMtz+7H9UDKvnBn9/i/Q3bog5HJC9F0Y030WIOSR8/mVkZ8CDQAzgTGA98j6Dx/NYGT+Y+1d2r3b26e/fu6UUsOaekuIgbTxuNu/O9h/7NXq0hIpJx2U4gGwhqIfEqSVwziXUuMAE43t3vc/eX3P2XwHeAi8zswIxGKjmvX5dyfnTicGYtWsc9s5ZEHY5I3sl2AqkhaAeJNxx4O0nZA4AN7v5e3PZ/hO/7NzM2yUOnHdyPo4b14OdPzuO9NVuiDkckr2Q7gTwOjDWzQfUbzKyKoIvu40nKrgIqzWxw3PZPh+9qLZVPCGbuPYB2bYr5jpbDFcmobCeQO4AlwDQzm2xmk4BpwHLg9vqDzGyAmdWZ2ZSYsncTNJz/zczONrPPmNn3gF8Ccwm6Aot8Qo+OZVw9aQSvL9/IvXqUJZIxWU0gYVfbo4AFwL3A/cBi4Ch3j32+YEBxbHzuvgQYC7xOMHr9bwQDE6cCn3V3zaInDZp0YB+OGNKN66fPZ+Wm7VGHI5IXzD21Kr2ZdQSOB/oDZXG73d1/kqHYMq66utrnzJkTdRgSkWXrtvG5m2Zw5JDuTD2rOupwRHKCmc1194T/w6Q0F5aZHQY8AXRu4BAHWm0CkcLWv2s5lx+zHz9/ch5PvbWKiSN7RR2SSE5L9RHWTQRtGAcDZe5eFPcqznSAIpl07uED2b93R656vIatO+uiDkckp6WaQPYHfujuc91dy79JziktLuLak0eyavMObnsxvke4iKQi1QSyDGjbEoGIZMuYAZWcPLoPU19exLJ1muZEJF2pJpCrgSvChnSRnHXFcftTUmT87G9awVAkXakuKHUi0BNYbGazCCZBjOXufnZGIhNpQb06lXHJZwZz/fT5vLpwLYcO7hZ1SCI5J9UayOEEPa02E0xJckSCl0hOOPfwgfTr0o6rn3hbI9RF0pBSAnH3gUleg5J/ikjrUFZazJXH7c/81bU88prWURdJVRTTuYu0GseN7MWB+3TipmcWsGO3lsAVSUXKCcTMys3sUjN7yMyeM7M/mdnXtSa55CIz4/sTh/HBph3cN3tp1OGI5JSUEoiZ9QJeA34NVAPlBIMKbwHmmlnPjEco0sIOHdyNI4Z045YXFrJ5x+6owxHJGanWQP6HYPGnI8I2j3HuPpCgcb0z8IsMxyeSFd+fOIyN23Zzx0uLog5FJGekmkCOA6509/+YOt3dXwV+CJyQqcBEsmlk306cOKo3v3t5Meu27Iw6HJGckGoCaQ980MC+98P9Ijnp8mP2Y0fdHn73yuKoQxHJCakmkPnAmQ3s+wowr3nhiERncI/2nHBAb/7w6hI2bNVUbyLJpJpAfgl8ycyeNbOvmdlxZnaOmU0HzgCuz3yIItnzjaOGsHXXHn4/U7UQkWRSmsrE3e8Lu+teA/wuZtdq4CJ3/2MmgxPJtqG9OnDcyF7cPXMJ5x0+iE7lpVGHJNJqpTwOxN2nAn34eCqTEUBfd78jw7GJROLSowZTu7OOu15VLUSkMWmNRHf3ve7+jrvPDN+1HrnkjRF9OnHM/j24+9UlbN+l0ekiDUn6CMvMzgL+6u7rwp8b5e5/yEhkIhG6cPy+PPvbWTw8dzlnjquKOhyRVqkpbSB3A2OBdeHPjXFACURyXvWASkb368zvXlnMGZ8eQHGRRR2SSKvTlEdYA4HXY35u7KXZeCUvmBnnHzGIpeu28czbq6MOR6RVSloDcfeliX4WyXfHjuhJvy7tuOPlRUwc2SvqcERanVQnU9xjZoc0sG+MmanFUfJGSXER5x42kLlLNzB36YaowxFpdVLthdXYg+BigjYQkbxxSnU/OpaVcPerS6IORaTVaVICMbMiMyuuLxP+HvuqIJhocW2LRSoSgYq2JZxS3Y+n3lrJh7U7og5HpFVJmkDM7MfAbmAXQQ1jZvh77GszMAV4qAmf18/MHjazTWa22cweNbP+TQ3YzPYPF7Naa2bbzWy+mV3W1PIiqfryp/uze4/zp38ujzoUkValKd14XwzfjSBJ3Ekw826sncDbwF8a+6BwGpTnw+PPJkhI1wIvmNkod9+apHx1WP5F4DxgEzAEzQIsLWhQ9/YcMaQbf/z7Mi4avy8lxVoJWgSa1gtrBjADwMwcuMPdG5rSPZnzCbr6DnX3heFn/ht4F7gQuKGhgmZWBNwDPOfun4/Z9UKasYg02VfGDuDCe+fy/LwP+dwI9cgSgRQb0d396mYkD4BJwOz65BF+5mKCx2KTk5SdAAynkSQj0lKOHtaD3p3KuFfrpot8JKXZeAHMrAfwJWAoUBa329393EaKjwCmJdheA5yS5NSHh+9lZjYbGANsAB4Avu/u25PFLpKukuIivnRIf254ZgFL1m6lqltF1CGJRC6lBGJmQ4HZBF12Kwh6XXUJf99A0CbRmC7hcfHWE6y13pg+4fuDwC3AFUA1wdTy/YDPJypkZhcAFwD079/ktnqRTzjt4H7c9OwCHpq7nO8dOyzqcEQil2pr4PXAP4CeBI3qxwHtCBq0t9HAP+JxEo0VacpEQ/Wx3ufuU9z9RXf/JXA1cLKZDU94Mvep7l7t7tXdu3dvwmlEEuvZsYzx+3Xnkbkr2LNXQ55EUk0gBwP/S9CLCqDI3evc/ffAb4CbkpTfQFALiVdJ4ppJrHXh+zNx258O30cnKS/SbKdU92PV5h28/O6aqEMRiVyqCaQ9sD5c/2MT0C1m3xyCBNOYGoJ2kHjDCboBJysLn6zB1NdetCaJtLij9+9B5/JSHpob35NdpPCkmkCWAPV9GOfznw3fJwIbk5R/HBhrZh/N2mtmVcBh4b7GPElQ85kYt/3Y8H1OkvIizda2pJiTR/flmZrVbNy2K+pwRCKVagJ5Bvhs+PMNwDnhSPAa4DLg90nK30GQhKaZ2WQzm0TQK2s5cHv9QWY2wMzqzGxK/TZ3XwdcB1xkZj8zs2PM7AqCwY33xHYNFmlJXxyzD7v27OXxN5rTo10k96XajfdKoC2Au//JzLYDpwHlwM0ECaJB7r7VzI4CbgTuJXj89BxwubtviTnUCHp2xSe4a4Ba4OvAd4GVBA37P0nxe4ikbWTfTgzr1YFHX1vBWVqtUApYSgnE3XfycQM67v4E8ESKn7EM+EKSY5aQoGeWuztBzUeDCSVSk0f35RdPzWP5+m3061IedTgikUh1PZApZvbFBvb1jX3kJJLPTjqwN4AeY0lBS7UN5CrgQTO7PsG+fYAfNzsikRywT2U5YwZU8oQSiBSwdKYVvRW4NJyGvV2mAxLJFSeN6s28VbUsWF0bdSgikUgngdwHHEPQ9fYlM9PUpFKQjh/VmyJDtRApWGktbODuM4FPE0xj8k8zG5XRqERyQI8OZYzbtytPvPEBQf8OkcKS9so4YU+pccCbwCvASRmKSSRnnDSqD0vWbaPmg81RhyKSdc1aWs3dawlGoN8F/HdGIhLJIZ8d3pMig6ffXh11KCJZl2oCOQd4L3aDu+9198uAMwgG+okUjK7t21I9oAtP16yKOhSRrEt1RcJ7wilFEu17wN2vzkxYIrnjcyN6Mm9VLcvWbYs6FJGsSjoS3czOAv7q7uvCnxvj7n5vZkITyQ2fG96La//6Dk+/vYrzjhiUvIBInmjKVCZ3A2MJ1uO4O8mxTjDHlUjB6N+1nGG9OvB0zWolECkoTUkgAwkmLaz/WUTifG5EL255/l3WbtlJt/Ztow5HJCuStoG4+1J33xXzc6Ovlg9ZpPX53PCe7HV47h31xpLCkepkinvM7JAG9o0xsz2ZCUskt4zo05E+ncp4ft6HUYcikjWpduP9xBTrMYr55HKzIgXBzBg/tAczF65j9x6triyFoUkJxMyKzKy4vkz4e+yrAjgOWNtikYq0cuP3686WnXW8tnRD1KGIZEXSBGJmPwZ2A7sIahgzw99jX5sJlpZ9qMUiFWnlDh3clZIi46V310QdikhWNKUX1ovhuxEkiTuB9+OO2Qm8DfwlY5GJ5JiOZaUcNKCSGQvW8L1jh0UdjkiLS5pA3H0GMAPAzBy4w901f7VIAuP368710+ezpnYn3TuoO6/kt1SnMrna3T8I2z1Gmtn4sP1DRAgSCMDLeowlBSDl2XjN7BJgFfBv4HlgaLj9MTP7ZmbDE8ktw3t3pFv7NsxYoAQi+S/VcSDnAzcDjwGn8p/del8GvpCxyERyUFGRcei+3Xj1vXVaZEryXqo1kG8Dv3L3C4A/x+2bR1gbESlk4/btypranby3ZmvUoYi0qFQTyEBgegP7tgKdmxWNSB4YN6grALMXJVz5QCRvpJpA1gJVDewbCqxoVjQieWBA13J6dypjlhKI5LlUE8gTwBQzi52z2s2sG/AtgrYRkYJmZowd1JW/L1I7iOS3VBPIDwkGDb4FPEswMv3XwDvAHrSkrQgQPMZau2UX7364JepQRFpMquNA1gHVwHVAKcH66CXALcA4d9+U7DPMrJ+ZPWxmm8xss5k9amb9Uw3czK40MzezV1ItK9LSxu0btIPMek+PsSR/pTwOxN1r3f0n7n64u+/n7uPCAYabk5U1s3KCsSPDgLOBM4EhwAupDEgMH6H9ANDc2dIq9etSTp9OZczRxIqSx5oyF1YmnQ8MAoa6+0IAM/s38C5wIXBDEz/nNuB+gob7bH8HkSY5aEClZuaVvJbOSPSzzewpM3vbzBbFvd5LUnwSMLs+eQC4+2KCGX4nN/H8ZwAHAVemGrtINo0ZUMmKjdtZuWl71KGItIiU/no3sx8BVxM0or9O0KCeihHAtATba4BTmnD+SuBG4L/cfb1ZY+tbiUTroP6VALy2dCMnjGoXcTQimZfq459zgZvd/Vtpnq8LkKhOvx6obEL564EFwN1NPaGZXQBcANC/f8pt9SJpG96nI2WlRcxduoETRvWOOhyRjEv1EVZXgrEgzZGoY3zSqoSZHQGcBVzsKXSud/ep7l7t7tXdu3dPIUyR5iktLmLUPp2Zu0ztIJKfUk0gM4ADm3G+DQS1kHiVJK6ZxLqdcDErM+tsZp0JalDF4e9afEFanTEDKqlZsYkdu/dEHYpIxqWaQC4HzjGzs8ysW4K10ZN9Xg1BO0i84QQrGjZmf+AigkRT/zoMGBv+fHEK30MkK8b0r6Rur/PG8o1RhyKScam2gSwI3+9qYL8n+czHgV+a2SB3XwRgZlUEieCKJOf+TIJtNwHFwDeAhQn2i0RqdP/OAPz7/U18OpxkUSRfpJpAriFxG0ZT3QFcCkwzsx+Gn/UTYDnBIyoAzGwAwSj3a9z9GgB3fzH+w8xsI1CSaJ9Ia9CtfVv6dCrjzRVJJ2kQyTkpJRB3v6o5J3P3rWZ2FEFX3HsJGs+fAy5399hJg4ygZpHyOBWR1mZk3068pQQieSjro7jdfRlJVi509yU0oWeWu0/ITFQiLeeAvp14+u3V1O7YTYey0qjDEcmYVAcSTmlk915gE/Cau89sVlQieWTkPp0AqPlgM2PVDiJ5JNUayFUE7RaJagf1293MZgEnNGV2XpF8d0DfIIG8tWKTEojklVTbGPYn6O30HWAAUBa+fy/cfihwenjczzIXpkju6ta+Lb3VkC55KNUayK3A79z9xphty4FfmVkx8FN3P9rMBhJ0rb0kQ3GK5LSRfTspgUjeSbUGMg54rYF9rxEM6gOYA/RINyiRfLN/744sWbtVI9Ilr6SaQDYBRzew75hwPwSPtpIuMCVSKIb16sBeh4Va4lbySKqPsH4PXGlmHYCHCVYE7EEwFftFBEvdAnyaYMp3EQH269kBgPmrahkZNqqL5LpUE0h9N97LgK+HPxuwlSB51O//K/Bgs6MTyRNVXctpU1LEgtW1UYcikjGpjkTfC/zQzH4JHAD0BlYCb7r7xpjj/pHJIEVyXUlxEYO7t2feKiUQyR9pjUQPk8XLmQ1FJL8N7dWB2YvWRR2GSMaklUDCpWWHEDSW/wd3f6m5QYnko6G9OvDnf61g07bddCrXlCaS+1KdyqSMoCH9VBqeq6q4uUGJ5KOhYUP6gg9rObgq0bpqIrkl1W68PwImAGcTJJBLgfOAVwimXz8xk8GJ5JP9en3cE0skH6SaQL5AsCbIA+Hvf3f3u9x9PPAGMDGTwYnkk94dyygrLWLJ2q1RhyKSEakmkP5AjbvvAXYDFTH7fg+clqnARPJNUZFR1bWCxUogkidSTSDrgPbhz8uBA2P2dQPaZSIokXw1qLsSiOSPVHthzQY+BTwJPAL8JByVXkcwQ+8rmQ1PJL8M7FbB0zWr2b1nL6XFWnBTcluqCeQXBI+xAK4FBhO0iRQTJJeLMxeaSP6p6lpB3V7n/Q3bGditInkBkVYs1ZHocwhm2sXda4EvmFlboK27a/JEkSQGdQ+SxuK1W5RAJOclTSBmdlRTPsgsGBbi7s83MyaRvDWwW9CEuGjNVo4aFnEwIs3UlBrIswTL1ULDgwc/Ws4WDSQUaVBleSmd2pWqIV3yQlMfYdUSNJo/QjDzroikwczo16Ud72/YHnUoIs3WlATyGeAsgkGEpwB/Bu7RoyqR9PTt3I731ujvMMl9SfsRuvsMdz8X6EWwaFQPYLqZLTOz68xs/5YOUiSf9O1czooN23H35AeLtGJN7oju7jvc/Y/ufhxBV96bgeOBt8zslpYKUCTf9K1sx/bde9iwbXfUoYg0S7ojmdYBS8KXA5UZikck7/XtHEzYsELtIJLjUkogZnaYmf2WYBXCe4AtwAnAmSl8Rj8ze9jMNpnZZjN71Mz6N6FctZlNNbN5ZrYtfIR2v5kNTOU7iERtn8owgWzcFnEkIs3TlHEggwkSxFeAKuAl4LvAQ+6+JZWTmVk58Dywk2BKeCcY0f6CmY1y98ZaFk8HRgC/BmqAvgTTy88xs9HuvjyVWESiUp9A1BNLcl1TemEtADYDjxKs/bE03N7DzHrEH+zuixr5rPOBQcBQd18IYGb/Bt4FLgRuaKTsL9x9TewGM5sJLA4/d0oTvotI5Dq1K6WiTbESiOS8po4D6Qh8laDWkExjAwknAbPrkweAuy8OE8FkGkkg8ckj3LbUzNYQ1EZEcoKZ0beyHSs2KoFIbmtKAjkng+cbAUxLsL2GYIxJSsIuxD2Ad5oZl0hW9e3cTo3okvOSJhB3vyeD5+sCbEiwfT0p9uQysxLgt8Aa4M5GjrsAuACgf/+kbfUiWdGrUxlvrtD8o5LboliQINHoqYbm2GrMLcChwFfcPVFSCk7mPtXdq929unv37mmcRiTzurdvy/qtO9mzV4MJJXdlO4FsIKiFxKskcc0kITO7jqBW8TV3fzpDsYlkTfcObdnrsG7rzqhDEUlbthNIDUE7SLzhwNtN+QAz+wFwBXCZu9+bwdhEsqZ7h7YArK3dFXEkIunLdgJ5HBhrZoPqN5hZFXBYuK9RZvZNgnEjP3D337RUkCItrT6BrNmiGojkrmwnkDsIpj+ZZmaTzWwSQa+s5cDt9QeZ2QAzqzOzKTHbTgduAp4CnjezsTGv4dn8EiLN1a19mEBqlUAkd6W6JnqzuPvWcIXDG4F7CRrPnwMujxvVbgTjSWIT3MRw+8TwFWsGMKGFwhbJOCUQyQdZTSAA7r6MYG2Rxo5ZQlzPLHf/KsFgRpGcV9G2hIo2xUogktOi6MYrIgTtIGoDkVymBCISkW7t27KmdkfUYYikTQlEJCKdy9uwaXtd1GGIpE0JRCQindqVsmmbxoFI7lICEYlI5/JSNm3XsraSu5RARCLSuV0pW3ftYfeevVGHIpIWJRCRiHQqLwVQLURylhKISEQ6tQsSyMZtSiCSm5RARCJSn0BUA5FcpQQiEpHO5W0A2LRdPbEkNymBiESkY1kwk5BqIJKrlEBEItK+bZBAtu7cE3EkIulRAhGJSEWYQLbt0mh0yU1KICIRaVdaDMAW1UAkRymBiESkqMioaFPMtp2qgUhuUgIRiVB52xK26hGW5CglEJEIVbQp1iMsyVlKICIRalNSxO46zYUluUkJRCRCJUVF1O1VApHcpAQiEqHSYmP3Ho86DJG0KIGIRKikWDUQyV1KICIRKilSDURylxKISITalBRpQSnJWUogIhEqKTLqVAORHKUEIhKhkmLVQCR3KYGIRKi02KjbqxqI5CYlEJEIlRQVUacaiOSorCcQM+tnZg+b2SYz22xmj5pZ/yaWLTOz681spZltN7NZZnZkS8cs0lJKNA5EclhWE4iZlQPPA8OAs4EzgSHAC2ZW0YSPuBM4H5gCnAisBKab2egWCVikhZVqJLrksJIsn+98YBAw1N0XApjZv4F3gQuBGxoqaGYHAmcAX3P3u8JtM4Aa4BpgUsuGLpJ5pSXqhSW5K9uPsCYBs+uTB4C7LwZmApObUHY38GBM2TrgAeBYM2ub+XBFWlZJURG71AYiOSrbNZARwLQE22uAU5pQdrG7b0tQtg0wOPxZJGeUFhtbdtbx2RtmRB2K5LHK8jb86aJxGf/cbCeQLsCGBNvXA5XNKFu//xPM7ALgAoD+/ZvUVi+SNSeO6sMHm3bgrsdY0nI6lpW2yOdmO4EAJPo/xZpQztIp6+5TgakA1dXV+r9UWpUD+3Xm1jMOijoMkbRkuw1kA4lrCpUkrl3EWt9I2fr9IiKSJdlOIDUEbRnxhgNvN6HswLArcHzZXcDCTxYREZGWku0E8jgw1swG1W8wsyrgsHBfsrKlxDS2m1kJcBrwtLvvzHi0IiLSoGwnkDuAJcA0M5tsZpMIemUtB26vP8jMBphZnZlNqd/m7q8TdOG9yczOM7OjCbrwDgR+nL2vICIikOUE4u5bgaOABcC9wP3AYuAod98Sc6gBxQniOwe4C7gW+CvQD5jo7q+1cOgiIhIn672w3H0Z8IUkxywhQe8qd98OfDt8iYhIhDQbr4iIpEUJRERE0mKFNALWzNYAS1Ms1g1Y2wLh5DNds9ToeqVO1yw1zbleA9y9e6IdBZVA0mFmc9y9Ouo4comuWWp0vVKna5aalrpeeoQlIiJpUQIREZG0KIEkNzXqAHKQrllqdL1Sp2uWmha5XmoDERGRtKgGIiIiaVECERGRtBRsAjGzfmb2sJltMrPNZvaomTVpyUIzKzOz681spZltN7NZZnZkS8ccpWZeL2/gNbqFw46Mme1jZr8J741t4fetamLZgru/oNnXrBDvsS+a2SNmtjS8T+ab2XVm1qEJZTNyjxVkAgnXFHkeGAacDZwJDAFeMLOKJnzEncD5wBTgRGAlMD1fb9YMXC+Au4Fxca8FGQ+29RgMnEqwUNrLKZYtqPsrRnOuGRTePfZdYA/w38BE4DbgYuAZM0v2b3tm7jF3L7gXcFl44QfHbBsI1AHfTlL2QIKldc+J2VYCzAcej/q7tbbrFR7rwLVRf48sX7OimJ/PC69BVRPKFdz91dxrFh5fiPdY9wTbzgqvxVGNlMvYPVaQNRBgEjDb3T9axdDdFwMzgclNKLubYG2S+rJ1BGuTHGtmbTMfbuSac70KkrvvTbNoId5fQLOuWUFy9zUJNv8zfO/bSNGM3WOFmkBGAG8l2F5DsERusrKL3X1bgrJtCKrh+aY516vexWa2M3y2/byZHZG58PJKId5fmaJ7DMaH7+80ckzG7rFCTSBdCJ6zxlsPVDajbP3+fNOc6wVwH/B14BjgAqAr8LyZTchQfPmkEO+vTCj4e8zM+gLXAM+6+5xGDs3YPZb1BaVakUQjKD+xiFUDx6RbNpel/Z3d/cyYX182s2kENZprgcMzEFs+KdT7q1kK/R4zs/YEy4PXEazc2ujhZOgeK9QayAYSZ9lKEmfmWOsbKVu/P98053p9grvXEixJfHAz48pHhXh/ZVwh3WNmVgY8DgwCjnX395MUydg9VqgJpIbgOWC84cDbTSg7MOzaGl92F7Dwk0VyXnOuV0Ma+iuo0BXi/dVS8v4eM7NS4BHgEOB4d3+zCcUydo8VagJ5HBhrZoPqN4QDlg4L9yUrWwqcElO2BDgNeNrdd2Y82ug153p9gpl1BE4A/p6pAPNIId5fGVcI91g41uN+4GhgsrvPbmLRzN1jUfdljqj/dAVBln2ToBvqJOANYBHQPua4AQTPFKfElX+A4NHNeeF/vIeBHcBBUX+31na9CAY73QGcAUwgGIj4JsFfOkdE/d1a+Lp9MXzdRvCX8MXh7+N1f2XumhXqPRZzja4Fxsa99snGPRb5RYjw4vcnqPptBmqBx4gbtARUhf+Brorb3g64AVgVXvS/AxOi/k6t8XoBJxGMF1lL0Pd8HcFfQIdE/Z2ycM28gdeLur8yd80K9R4DljRyva7Kxj2m6dxFRCQthdoGIiIizaQEIiIiaVECERGRtCiBiIhIWpRAREQkLUogIiKSFiUQyVmNLGMa+1oSdZzZZmZVZnZV7MwBIi2hkGfjldw3Lu73PxOMkL8qZlshTv1RBfwYeIVgtgCRFqEEIjnL4+b+MbOdwNr47fnAzNp6xPNgtYYYpHXRIyzJa2Y20MzuN7M14Wp1r5vZ5+OOuSp83DXMzKab2VYzW2Zm54T7zzSzeWa2xcxeMLN948ovMbP7zOx8M1toZjvM7DUz+0yCeMab2XNmVhueZ7qZjYw75kUze8XMTjKzf4WJ8evhvkvNbJaZrTezjWY228xOiCk7AXgh/PWZmEd5E8L9bmZXxZ2vKtz+1Zhtd5vZ+2Y2zsxeNbPtwP+E+7qZ2W1mtiK8pvPM7IKm/1eRfKEEInnLzPoRzPFzIPAtgkkgXwMeMbNJCYo8RLCGxMnAXOD3ZvYzggn9riBYqGco8McEZccD3wZ+AJxO8OjsSTMbGhPPCcBzwBbgKwST/3UgWACpX9zn7Qf8GvgNcGxYDoLHU78jmEn1NGAO8BczOy7c/xpwSfjzNwke840Lt6eqE8Gke/8HHAf8MZzldibBTLdXhe9PALeZ2TfSOIfksqgnBNNLr0y9CCaXuy/m9zuBNUDXuOOeAV6P+f0qggnnzorZVkkwi+k6oGPM9m+Gxw6IO+8uoH/Mtg4EC/PcG7NtIfBcXCwdCSYBvClm24vAXmB0ku9bRPAY+mlgWsz2CWGMxyQok2hivapw+1djtt0dbpscd+yPCCbfGxK3/Y7we5REfR/olb2XaiCSzyYCfwM2mVlJ/QuYDhwY/jUd68n6H9x9A/AhMNvdN8ccMy98j68xzHb3ZTHl61fEGwdgZkOAfYH742LZBswCjoz7vCXu/nr8FzKzMWb2FzNbTZDgdgOfJagZZVod8Je4bRMJanWLE1zTrgSLEkmBUCO65LMewFnhK5GuBNPT14tfnndXA9sAyuK2r07w+auBvjGxQFArujPBscvifl8Zf0D4mOs5glUgvxGWqQN+Auyf4DOb60N33xO3rQcwmCBxJdK1BeKQVkoJRPLZOuBl4BcN7P8gg+fq2cC2FTGxAFwJPJvg2F1xvydaZ2EiQbvEqR6z7nWCpUkbsxNoE7etoX/0E8WwjqBmdlkDZeanEIvkOCUQyWdPETxCqnH37S18rrFm1s/dlwOYWQeCBua/hvvnE7SVjHD3n6d5jvpE8dFf/2a2H8HSwu/HHFff1bZdgs9YCoyM23ZCguMa8hRh7cfdP0yhnOQhJRDJZ1OAfwAvmdktBP+AVxL8AzrI3b+WwXOtBp4Ou8juBL5PsBTwTwDc3c3sEmCambUB/kTQ6NwTOJTgH+QbkpzjWYJHVn8ws18BvYGrCR5lxbZnLgiP+5qZrQ/jmR+2yzwA/NDMfgDMBo4AvpTC97yRoPfXy2Z2I0FirACGESwfOzmFz5Icp0Z0yVtho3Y1wej0nxH0vrqNoMvt8xk+3QzgV+F5HiRoIznO3RfExPM3gsbyCoKuuNMJxlb0ImhIb5S71wBfJljn+nHgvwi6F78Ud9w64FKC7sszgH8CY8Ld1wG3hPsfI2g7ObOpX9LdNxEkvL8RJMnpwO+ByXw8/kQKhJa0FWmmcL6tV9z9K1HHIpJNqoGIiEhalEBERCQteoQlIiJpUQ1ERETSogQiIiJpUQIREZG0KIGIiEhalEBERCQt/x/Q6ClmW7iDFQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "from math import tanh,cosh\n", "from numpy import linspace\n", "import matplotlib.pyplot as plt \n", "\n", "# Constants\n", "Tmax = 2.0\n", "points = 1000\n", "accuracy = 1e-6\n", "\n", "# Set up lists for plotting\n", "y = []\n", "temp = linspace(0.01,Tmax,points)\n", "\n", "# Temperature loop\n", "for T in temp:\n", " m1 = 1.0\n", " error = 1.0\n", "\n", " # Loop until error is small enough\n", " while error>accuracy:\n", " m1,m2 = tanh(m1/T),m1\n", " error = abs((m1-m2)/(1-T*cosh(m2/T)**2))\n", " y.append(m1)\n", "\n", "# Make the graph\n", "plt.rc('font',size=16) # set the font size\n", "plt.plot(temp,y)\n", "plt.xlabel('Temperature')\n", "plt.ylabel('Magnetization')\n", "\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercise 2 - Newton's method\n", "\n", "1. We can now return to our earlier example of the root of the function $f(x)=x - e^{1-x^2}$, for which the simple relaxation method failed. For the Newton method we need to add the derivative, which we have already derived earlier $f'(x)=1+2xe^{1-x^2}$. Modifying your relaxation method code to incorporate the expression for the Newton method.\n", "\n", "Use Newton's method to calculate the inverse (or arc) hyperbolic tangent of a number $u$. By definition, $\\tanh^{-1}{u}$ is the number $x$ such that $u=\\tanh{(x)}$. To put that another way, $x$ is a root of the equation $\\tanh{(x)}-u=0$. Recalling that the derivative of $\\tanh{(x)}$ is $1/\\cosh^2{(x)}$ the next Guess in Newton's method becomes $$x'=x-(\\tanh{(x)}-u)\\cosh^2{(x)}.$$ \n", "\n", "2. Write a function that calculates $\\tanh^{-1}{(u)}$.\n", "3. Plot $\\tanh^{-1}{(u)}$ from -1 to 1.\n", "\n", "\n", "## Talking points\n", "\n", "1. What do you observe?\n", "2. How quickly does Newton's method find the right solution?\n", "3. Does your function $\\tanh^{-1}{(u)}$ give the right solution?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercise 3 - Golden ratio search\n", "\n", "The Buckingham potential is an approximate representation of the potential energy between atoms in a solid or gas as a function of the distance $r$ between them: $$V(r)=V_0\\left[\\left(\\frac{\\sigma}{r}\\right)^6-e^{-r/\\sigma} \\right] .$$ \n", "\n", "1. Plot the Buckingham potential for $\\sigma$=1.\n", "2. Complete the golden ratio example program below to find the minimum of the Buckingham potential.\n", "3. Check your computational against the analytic solution.\n", "\n", "## Talking points\n", "\n", "1. What do you observe?\n", "2. What can you say about the Buckingham potential?\n", "3. How does the number of iterations depend on the specified accuracy?\n", "\n", "Example program to be completed:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "from math import exp,sqrt\n", "\n", "# Constants\n", "sigma = 1.0 # Value of sigma in nm\n", "accuracy = 1e-6 # Required accuracy in nm\n", "z = (1+sqrt(5))/2 # Golden ratio\n", "\n", "# Function to calculate the Buckingham potential\n", "def f(r):\n", " return (sigma/r)**6 - exp(-r/sigma)\n", "\n", "# Initial positions of the four points\n", "x1 = sigma/10\n", "x4 = sigma*10\n", "x2 = x4 - (x4-x1)/z\n", "x3 = x1 + (x4-x1)/z\n", "\n", "# Initial values of the function at the four points\n", "f1 = f(x1)\n", "f2 = f(x2)\n", "f3 = f(x3)\n", "f4 = f(x4)\n", "\n", "# Main loop of the search process\n", "\n", "# ADD THE GOLDEN RATIO SEARCH PROCEDURE HERE. YOU WILL NEED AN IF STATEMENT AND A STOPPING CONDITION.\n", "\n", "\n", "# Print the result\n", "print(\"The minimum falls at\",0.5*(x1+x4),\"nm\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercise 4 - Gauss-Newton method and gradient descent\n", "\n", "Return to the Buckingham potential from the previous exercise: $$V(r)=V_0\\left[\\left(\\frac{\\sigma}{r}\\right)^6-e^{-r/\\sigma} \\right] .$$ Find the minimum of the Buckingham potential for $\\sigma$=1:\n", "1. For the Gauss-Newton method. Start from $r=\\sigma$.\n", "2. For gradient descent.\n", "3. For gradient descent with numeric 1st derivative.\n", "\n", "## Talking points\n", "\n", "1. What do you observe?\n", "2. What happens when you start from $r=4\\sigma$ and why?\n", "3. What is a good value for $\\gamma$ in gradient descent?\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.9" } }, "nbformat": 4, "nbformat_minor": 4 }