{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Solve the first order differential equations\n", "$\\frac{dy}{dt} = f(t,y)$, $y(t_0) = y_0$ and $a \\le t \\le b$.\n", "\n", "For example,\n", "$\\frac{dy}{dt} = -5 t y^2+\\frac{5}{t}-\\frac{1}{t^2}$, $y(0) = 1$ and $1 \\le t \\le 25$" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "collapsed": true }, "outputs": [], "source": [ "%plot gnuplot" ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "collapsed": true }, "outputs": [], "source": [ "function [y]=forwardeuler(a, b,h, y0, f)\n", "% function [y]=forwardeuler(a, b,h, y0, f)\n", "% Forward Euler Method\n", "% end points a and b\n", "% time step h\t\n", "% initial condition y0\n", "% input function f\n", "\n", "y(1)=y0;\n", "\n", "t=a : h : b;\n", "n=length(t);\n", "\n", "for i=1:n-1\n", " y(i+1)=y(i)+h*f(t(i),y(i));\n", "end\n", "end" ] }, { "cell_type": "code", "execution_count": 23, "metadata": { "collapsed": true }, "outputs": [], "source": [ "function [ y ] = rk4( a,b,h,y0, f)\n", "% function [ y ] = rk4( a,b,h,y0, f)\n", "% Runge-Kutta fourth order method\n", "% end points a and b\n", "% time step h\t\n", "% initial condition y0\n", "% input function f\n", "\n", "y(1)=y0;\n", "\n", "t=a : h : b;\n", "n=length(t);\n", "\n", "for i=1:n-1\n", " f1=f(t(i),y(i));\n", " Y2=y(i)+0.5*h*f1;\n", " f2=f(t(i)+h/2, Y2);\n", " Y3=y(i)+0.5*h*f2;\n", " f3=f(t(i)+h/2, Y3);\n", " Y4=y(i)+h*f3;\n", " f4=f(t(i)+h, Y4);\n", " y(i+1)=y(i)+(h/6)*(f1+2*f2+2*f3+f4);\n", "end\n", "\n", "end" ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "collapsed": false }, "outputs": [], "source": [ "clear all\n", "h= 0.1;\n", "a=1;\n", "b=25;\n", "y0=1;\n", "f=@(t,y) -5*t*y^2+(5/t)-(1/t^2);\n", "\n", "[y1] = forwardeuler(a, b, h, y0, f);\n", "[y2] = rk4(a,b,h,y0, f);\n", "\n", "t=a : h : b;\n", "t1=a : h*10 : b;\n", "yexact=@(t) 1./t;\n", "\n", "figure;\n", "\n", "subplot(1,2,1)\n", "plot(t,y1,t1,yexact(t1),'o')\n", "title('Forward Euler')\n", "\n", "subplot(1,2,2)\n", "plot(t,y2,t1,yexact(t1),'o')\n", "title('RK4')\n", "print -dpng ode.png" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Only Matlab" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\\begin{equation}\n", "\\begin{split}\n", "\\frac{dX}{dt} & = -X - Y,\\\\\n", "\\frac{dY}{dt} & = -pZ + rY + sZ^2 -Z^2Y,\\\\\n", "\\frac{dZ}{dt} & = -q(X+Z).\n", "\\end{split}\n", "\\end{equation}\n", "\n", "
\n", " \n", " | \n", "\n", " \n", " | \n", "
---|