{ "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": [ "\"Phase" ] }, { "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", "" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "function [x,y,z] = climate(p, r, s,q, initV, T, eps)\n", "% \n", "% INITV - initial point\n", "% T - time interval\n", "% EPS - ode solver precision\n", "%\n", "% Example.\n", "% [X Y Z]=climate(1,0.8,0.8,1.2)\n", "% plot3(X,Y,Z);\n", "\n", "if nargin<4\n", " error('MATLAB:lorenz:NotEnoughInputs','Not enough input arguments.'); \n", "end\n", "\n", "if nargin<5\n", " eps = 0.000001;\n", " T = [0 10];\n", " initV = [0.001 0.001 0.001];\n", "end\n", "\n", "if nargin<6\n", " eps = 0.000001;\n", " T = [0 100];\n", "end\n", "\n", "options = odeset('RelTol',eps,'AbsTol',[eps eps eps/10]);\n", "[T,X] = ode45(@(T,X) F(T, X, p, r, s,q), T, initV, options);\n", "Tscale = 10;\n", "Tdim = T*Tscale;\n", "figure(1)\n", "plot3(X(:,1),X(:,2),X(:,3),'r');\n", "axis equal;\n", "grid;\n", "title('Phase portraits');\n", "xlabel('X'); ylabel('Y'); zlabel('Z');\n", "hold on\n", "\n", "figure(2)\n", "subplot(3,1,1);\n", "plot(Tdim,X(:,1),'r');\n", "ylabel('X'); \n", "subplot(3,1,2);\n", "plot(Tdim,X(:,2),'g');\n", "ylabel('Y'); \n", "subplot(3,1,3);\n", "plot(Tdim,X(:,3),'b');\n", "xlabel('kyr'); ylabel('Z');\n", "\n", "x = X(:,1);\n", "y = X(:,2);\n", "z = X(:,3);\n", "return\n", "end\n", "\n", "function dx = F(T, X, p, r, s,q);\n", "\n", " dx = zeros(3,1);\n", " dx(1) = -X(1) - X(2);\n", " dx(2) = -p*X(3)+r*X(2) + X(3)*X(3)*(s-X(2));\n", " dx(3) = -q*(X(1)+X(3));\n", "end" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "[X Y Z]=climate(1,0.8,0.8,1.2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ " \n", " \n", " \n", " \n", " \n", " \n", "
\n", " \"Phase\n", "\n", " \"Climate\n", "
\n" ] } ], "metadata": { "kernelspec": { "display_name": "Octave", "language": "octave", "name": "octave" }, "language_info": { "file_extension": ".m", "help_links": [ { "text": "MetaKernel Magics", "url": "https://github.com/calysto/metakernel/blob/master/metakernel/magics/README.md" } ], "mimetype": "text/x-octave", "name": "octave", "version": "0.16.1" } }, "nbformat": 4, "nbformat_minor": 1 }