HW6_Problem2_demo.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import math\n",
    "import numpy.linalg as nla\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.integrate import solve_bvp, ode, odeint\n",
    "from scipy.optimize import least_squares\n",
    "plt.rc('text',usetex=True)\n",
    "from sympy import lambdify, simplify, Symbol, init_printing, log, solve, Function\n",
    "init_printing()"
   ]
  },
  {
   "cell_type": "code",
<<<<<<< HEAD
   "execution_count": 2,
=======
   "execution_count": null,
>>>>>>> facc6a54fef4bfd9faa0ac7d017a481958927f81
   "metadata": {},
   "outputs": [],
   "source": [
    "# Material properties:\n",
    "# mu = 1. \n",
    "# kapp = 20.\n",
    "# alph = 2.5\n",
    "\n",
    "mu=Symbol('mu')\n",
    "kapp=Symbol('kappa')\n",
    "alph = Symbol('alpha')\n",
    "lam1, lam2, lam3 = Symbol('lambda_1'), Symbol('lambda_2'), Symbol('lambda_3')\n",
    "lam = Symbol('lambda')\n",
    "R = Symbol('R')\n",
    "B = Symbol('B')\n",
    "f = Function('f', real=True)\n",
    "fp = Function('f_p', real=True)\n",
    "Jdet = lam1*lam2*lam3"
   ]
  },
  {
   "cell_type": "code",
<<<<<<< HEAD
   "execution_count": 3,
=======
   "execution_count": null,
>>>>>>> facc6a54fef4bfd9faa0ac7d017a481958927f81
   "metadata": {},
   "outputs": [],
   "source": [
    "# define the stored energy function in terms of the eigen-stretches\n",
    "\n",
    "W = mu/alph*(lam1**alph + lam2**alph + lam3**alph ) - mu*log(Jdet) + kapp/2*(Jdet - 1)**2"
   ]
  },
  {
   "cell_type": "code",
<<<<<<< HEAD
   "execution_count": 4,
=======
   "execution_count": null,
>>>>>>> facc6a54fef4bfd9faa0ac7d017a481958927f81
   "metadata": {},
   "outputs": [],
   "source": [
    "# compute the symbolic derivatives: \n",
    "    # useful in calculation of the eigen stresses\n",
    "    \n",
    "dW_dl1 = lambdify((lam1, lam2, lam3), simplify(W.diff(lam1)), 'numpy')\n",
    "dW_dl2 = lambdify((lam1, lam2, lam3), simplify(W.diff(lam2)), 'numpy')\n",
    "dW_dl3 = lambdify((lam1, lam2, lam3), simplify(W.diff(lam3)), 'numpy')\n",
    "\n",
    "# t1 = dW_dl1(R*fp+f, f, f)\n",
    "# t2 = dW_dl2(R*fp+f, f, f)\n",
    "# t3 = dW_dl3(R*fp+f, f, f)\n",
    "t1 = simplify(1./(lam2*lam3)*W.diff(lam1)).subs([(lam1,R*fp(R)+f(R)), (lam2, f(R)), (lam3, f(R))])\n",
    "t2 = simplify(1./(lam1*lam3)*W.diff(lam2)).subs([(lam1,R*fp(R)+f(R)), (lam2, f(R)), (lam3, f(R))])\n",
    "t3 = simplify(1./(lam1*lam2)*W.diff(lam3)).subs([(lam1,R*fp(R)+f(R)), (lam2, f(R)), (lam3, f(R))])"
<<<<<<< HEAD
=======
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t1_func = lambdify((lam1, lam2, lam3), simplify(1./(lam2*lam3)*W.diff(lam1)), 'numpy')  #takes three arguments (lambda_1, lambda_2, lambda_3)"
>>>>>>> facc6a54fef4bfd9faa0ac7d017a481958927f81
   ]
  },
  {
   "cell_type": "code",
<<<<<<< HEAD
   "execution_count": 5,
=======
   "execution_count": null,
>>>>>>> facc6a54fef4bfd9faa0ac7d017a481958927f81
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the second derivative f''(R), which will be helpful later on to impose BC and \n",
    "# solve the resulting 1st-order system\n",
    "\n",
    "fpp = simplify(solve(t1.diff(R)+2./R*(t1-t2),fp(R).diff(R))[0])\n",
    "fpp = simplify(fpp.subs(f(R).diff(R), fp(R)))\n",
    "fpp"
   ]
  },
  {
   "cell_type": "code",
<<<<<<< HEAD
   "execution_count": 8,
=======
   "execution_count": null,
>>>>>>> facc6a54fef4bfd9faa0ac7d017a481958927f81
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(-2.0*R**4*kappa*f(R)**4*f_p(R)**3*Derivative(f(R), R) - R**3*kappa*f(R)**5*f_p(R)**3 - 7.0*R**3*kappa*f(R)**5*f_p(R)**2*Derivative(f(R), R) - 2.0*R**2*kappa*f(R)**6*f_p(R)**2 - 8.0*R**2*kappa*f(R)**6*f_p(R)*Derivative(f(R), R) + 2.0*R**2*mu*(R*f_p(R) + f(R))**alpha*f_p(R)*Derivative(f(R), R) - 2.0*R**2*mu*f_p(R)*Derivative(f(R), R) - R*alpha*mu*(R*f_p(R) + f(R))**alpha*f(R)*f_p(R) - R*alpha*mu*(R*f_p(R) + f(R))**alpha*f(R)*Derivative(f(R), R) - R*kappa*f(R)**7*f_p(R) - 3.0*R*kappa*f(R)**7*Derivative(f(R), R) - R*mu*(R*f_p(R) + f(R))**alpha*f(R)*f_p(R) + 3.0*R*mu*(R*f_p(R) + f(R))**alpha*f(R)*Derivative(f(R), R) - R*mu*f(R)*f_p(R) - 3.0*R*mu*f(R)*Derivative(f(R), R) + 2.0*R*mu*f(R)**(alpha + 1)*f_p(R) - 2.0*mu*(R*f_p(R) + f(R))**alpha*f(R)**2 + 2.0*mu*f(R)**(alpha + 2))/(R**2*(R**2*kappa*f(R)**4*f_p(R)**2 + 2.0*R*kappa*f(R)**5*f_p(R) + alpha*mu*(R*f_p(R) + f(R))**alpha + kappa*f(R)**6 - mu*(R*f_p(R) + f(R))**alpha + mu)*f(R))\n"
     ]
    }
   ],
   "source": [
<<<<<<< HEAD
    "# Calculate the second derivative f''(R), which will be helpful later on to impose BC and \n",
    "# solve the resulting 1st-order system\n",
    "fpp = simplify(solve(t1.diff(R)+2/R*(t1-t2),fp(R).diff(R))[0])\n",
    "# fpp = simplify(fpp.subs(f(R).diff(R), fp(R)))\n",
    "print(fpp)"
=======
    "p = Symbol('p')\n",
    "# solve(p + )\n",
    "t1.subs([(f(R),lam), (alph,2.5)])"
>>>>>>> facc6a54fef4bfd9faa0ac7d017a481958927f81
   ]
  },
  {
   "cell_type": "code",
<<<<<<< HEAD
   "execution_count": 9,
=======
   "execution_count": null,
>>>>>>> facc6a54fef4bfd9faa0ac7d017a481958927f81
   "metadata": {},
   "outputs": [],
   "source": [
    "def ode_sys(r, z):\n",
    "    \"\"\"\n",
    "    Let's call the state variable as z = [f, f']\n",
    "        -- implies z' = [f', f''], where f'' comes from the ode\n",
    "    \"\"\"\n",
    "    mu, kappa, alpha = 1., 20., 2.5\n",
    "    f_r, fp_r = z\n",
    "    fpp_r = (-2.0*r**4*kappa*f_r**4*fp_r**4 - 8.0*r**3*kappa*f_r**5*fp_r**3 - 10.0*r**2*kappa*f_r**6*fp_r**2 + 2.0*r**2*mu*(r*fp_r + f_r)**alpha*fp_r**2 - 2.0*r**2*mu*fp_r**2 - 2*r*alpha*mu*(r*fp_r + f_r)**alpha*f_r*fp_r - 4.0*r*kappa*f_r**7*fp_r + 6.0*r*mu*(r*fp_r + f_r)**alpha*f_r*fp_r - 4.0*r*mu*f_r*fp_r - 2.0*r*mu*f_r**(alpha + 1)*fp_r + 2.0*mu*(r*fp_r + f_r)**alpha*f_r**2 - 2.0*mu*f_r**(alpha + 2))/(r**2*(r**2*kappa*f_r**4*fp_r**2 + 2.0*r*kappa*f_r**5*fp_r + alpha*mu*(r*fp_r + f_r)**alpha + kappa*f_r**6 - mu*(r*fp_r + f_r)**alpha + mu)*f_r)\n",
    "    return [fp_r, fpp_r]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bc_z(za, zb):\n",
    "    mu, kappa, alpha = 1., 20., 2.5\n",
    "    A = 0.95\n",
    "    B = 1.\n",
    "    f_A, fp_A = za\n",
    "    f_B, fp_B = zb\n",
    "    \n",
    "#     include the boundary condition in terms of the cauchy pressure\n",
    "    bc_inside = p + t1_func(A*fp_A + f_A, f_A, f_A)\n",
    "    bc_outside = t1_func(B*fp_B + f_B, f_B, f_B)\n",
    "    return [bc_inside, bc_outside]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = 0.95\n",
    "B = 1.\n",
    "pvals = np.linspace(0,0.09*10**6,100)\n",
    "rVals = np.linspace(A,B,1000)\n",
    "z_a = np.ones((2, rVals.size))\n",
    "z_a[1] = 2.\n",
    "rB = np.zeros(pvals.shape)\n",
    "rVals_plot = np.linspace(A,B,10**4)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for idx, p in enumerate(pvals):\n",
    "    z_a = np.ones((2, rVals.size))*(1+p**(1./4))\n",
    "    z_a[1] = 1./(1+p**(1./4))\n",
    "    res_a = solve_bvp(ode_sys,bc_z, rVals, z_a,tol=1.e-8, max_nodes = 1.e8,verbose=1)\n",
    "    rB[idx] = res_a.sol(rVals_plot)[0][-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.linspace(0,1,11)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "90000**(0.25)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "R = Symbol('R')\n",
    "A = Symbol('A')\n",
    "lam_z = Symbol('lambda_z')\n",
    "f = Function('f')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAJ8AAAA8BAMAAACKvXWTAAAAMFBMVEX///8AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEM3dMqtUiUSZIrt2Zu8ziBZcAAAACXBIWXMAAA7EAAAOxAGVKw4bAAADKElEQVRYCeWXP2gTURzHv5dees0lac52dEgcFAfRgIrrgTgJbQZFxaHBxUVIO2gRh2ao4B/UIKI4iAFHhwZEsaSFOEgopRIUHFwMRR3UIcU/rVWJ7727tu/dvcuf3k31DXe/3+f3/f1yuV6+7wq4VrjpYzVc44BBCfOFzvrqljQfkzBfKO2r293cV3AzXySc59sN9FZP8gA4d15QuOqiGujhgfIB/Uiuk+u/Lp5AXRGeDVpXjq9L3MFlHmk/oWOIEnYjYmmECkW1vg5IwOrJLEXydZfHk19J9pySEj2EC4imEM3TmAEa0HokxSLpIc3Tyk4iNilh/f1Z3DRxhykoUB5Ui6yurzImO6gZjsZwGxjTKGEDpxbmxqAZo2tg0Bw2rPo01yWGWpbLJ/EKoeY3StjAN9g+in3N9Bq4T+6vVR8xuDYh7KGVo2WyisC9cUsYKZe/lGeBU+htMLUNVvDSag4frliB+/iYQzFgyrTzEj1/glYnJ7YIUFYxu40myu54HdqhF1ZFPC5w6Q0gkbVzOlBfQjh9jQOf1ZUaTR/lyWfpcdMuCacdG1nkI5SHu+ycDpxfNtQ9FQ5cfVtlA58Br4F3dkU4KSkh3UjoQGG5gGZkBYGVxCsSSFHByV1g4UlN1LAHJloUoZ8slyHdzBoSm95UGtwFqGfqJAvQ//XYdzKwtf+3tCnu4qxQXSbntAsLoJVNCUKWLALt/L+VTbkHkt9tJO/GPGlhU7zMjnNZ0f8lEnjblESdzOCIBAvI26YEmZWE6rggwTxqYVO8zI7Jc9Pmj8xsStLpgdRlveRRsjGzKQ/fk3UuCv4vUTCb8vA9iRwjA4YMO5nU95wilucuSbEDyn3PIbLS5F8pdkCX7znqXBpa4pIgwlgqiCncDMXkkq0bDuyn6wA2v48IG9DWvVH/yTebbz4tE/8KcP0GhjMBzlP/AFHvt/3uP4m+qoSCHBhqkK/cZtPp6jITJegHu+poI741MzdtttF0VZ4oKqc72sQ6nbqX/JMx2qm4E90PgLz9BLf6yJaYKwQ3DxGyJU5U1MAmKkNk4PvKlcAGauPjNURmau0G/gPlH2VLrkkwzAAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$- \\frac{\\sqrt{\\frac{- A^{2} + R^{2} + a^{2} \\lambda_{z}}{\\lambda_{z}}}}{R}$$"
      ],
      "text/plain": [
       "      ____________________ \n",
       "     ╱    2    2    2      \n",
       "    ╱  - A  + R  + a ⋅λ_z  \n",
       "-  ╱   ──────────────────  \n",
       " ╲╱           λ_z          \n",
       "───────────────────────────\n",
       "             R             "
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "LHS = (1-lam_z*f(R)**2)/(1-lam_z*f(A)**2)\n",
    "RHS = A**2/R**2\n",
    "simplify(solve(simplify(LHS-RHS),f(R))[0]).subs(f(A)*A,Symbol('a'))"
   ]
  },
  {
   "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.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}