Skip to content

Instantly share code, notes, and snippets.

@marcosbarker

marcosbarker/lorenz_da_magic.ipynb

Created February 9, 2023 15:43
Show Gist options
  • Save marcosbarker/9fde9f3ccff6f00e36c30e5f0da592ee to your computer and use it in GitHub Desktop.
Save marcosbarker/9fde9f3ccff6f00e36c30e5f0da592ee to your computer and use it in GitHub Desktop.
Download ZIP
Lorenz_DA_Magic.ipynb
Display the source blob
Display the rendered blob
Raw
lorenz_da_magic.ipynb
{
"nbformat": 4,
"nbformat_minor": 0,
"metadata": {
"colab": {
"provenance": [],
"collapsed_sections": [
"BbCnpun9zGJE",
"Rnm9Hzi1M8jo",
"emSJfHDZFkll",
"mIYDfN2UA5Xh",
"OZ33hAOVBD-F",
"bGlX8KlFBn1k",
"CK6XUeR2Bwqk",
"Iniw4SxLB4pl",
"HFX7s8z-CBRG",
"e1lUNJMbCG87",
"Ixn-ri7YDJ1j",
"Me3ZlKoRDjiT",
"HKutT7p1ESh8",
"X2bPqh7MEh2c",
"apQJBpVywEew",
"-uKbAaNWwlKQ",
"m1-TAPgJ7INl",
"YHKwdGUn7qb-",
"8uGTh48s7aVz",
"iEyiPjj7eK9u",
"nyNs4Q-O2Qry",
"yASAkISp2a89",
"ltN_czPf2nXJ",
"txQqnLrk3gOx"
],
"include_colab_link": true
},
"kernelspec": {
"name": "python3",
"display_name": "Python 3"
},
"language_info": {
"name": "python"
}
},
"cells": [
{
"cell_type": "markdown",
"metadata": {
"id": "view-in-github",
"colab_type": "text"
},
"source": [
"<a href=\"https://colab.research.google.com/gist/marcosbarker/9fde9f3ccff6f00e36c30e5f0da592ee/lorenz_da_magic.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
]
},
{
"cell_type": "markdown",
"source": [
"# Program to perform data assimilation on the Lorenz equations\n",
"### (c) 2002 Data Assimilation Research Centre\n",
"### Original program by Matthew Martin"
],
"metadata": {
"id": "BbCnpun9zGJE"
}
},
{
"cell_type": "markdown",
"source": [
"# **Instala o Octave e o plugin para jupyter notebook**"
],
"metadata": {
"id": "Rnm9Hzi1M8jo"
}
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"id": "351NAv5Swr3S"
},
"outputs": [],
"source": [
"!apt install octave"
]
},
{
"cell_type": "code",
"source": [
"!pip install oct2py"
],
"metadata": {
"id": "Su0Mh-wEz4VM"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"source": [
"%load_ext oct2py.ipython"
],
"metadata": {
"id": "EbNdUds2E39h"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"# MAIN PROGRAM\n",
"\n"
],
"metadata": {
"id": "yE8qioOfLZnD"
}
},
{
"cell_type": "markdown",
"source": [
"## 0. PROGRAM's BRIEF DESCRIPTION"
],
"metadata": {
"id": "emSJfHDZFkll"
}
},
{
"cell_type": "code",
"source": [
"%%octave\n",
"\n",
"%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
"% Program to perform data assimilation on the\n",
"% Lorenz equations\n",
"%\n",
"%\n",
"% (c) 2002 Data Assimilation Research Centre\n",
"%\n",
"% Original program by Matthew Martin\n",
"%\n",
"% Change history:\n",
"% Changes to make code more robust and add further options\n",
"% (Amos Lawless)\n",
"% 23/01/03 Change to position windows correctly for any screen\n",
"% resolution (Amos Lawless)\n",
"% 16/07/04 Change to initial conditions and run time\n",
"%\n",
"% Structure of program:\n",
"% 1. Inputs\n",
"% 2. `True' solution\n",
"% 3. Background solution\n",
"% 4. Observations\n",
"% 5. Covariance matrices\n",
"% 6. Assimilation\n",
"% 7. Plot results.\n",
"%\n",
"% List of main variables\n",
"% h: Time step for Runge-Kutta scheme\n",
"% s: sigma coefficient in equations\n",
"% r: rho coefficient in equations\n",
"% b: beta coefficient in equations\n",
"%\n",
"% fmax: Total time of assimilation + forecast\n",
"% ob_f: Frequency of observations in time steps\n",
"% q: x-axis data values for plots\n",
"% sd: Variance of observation error\n",
"% tmax: Total time of assimilation\n",
"%\n",
"% R(3,3): Observation error covariance matrix\n",
"% B(3,3): Background error covariance matrix\n",
"% Qx(3,3,fmax): Model error covariance matrix for Kalman Filter\n",
"% (default constant in time)\n",
"%\n",
"% x(fmax,1), y(fmax,1), z(fmax,1):\n",
"% True state vectors of x,y,z\n",
"% xb(fmax,1), yb(fmax,1), zb(fmax,1):\n",
"% Background state vectors of x,y,z\n",
"% xob(tmax), yob(tmax), zob(tmax)\n",
"% Observations of x,y,z\n",
"%\n",
"% x_ob(3,tmax): [xob';yob';zob']\n",
"% x_sc(3,fmax): Analysis and forecast from SC scheme, with first\n",
"% index representing x, y or z variable\n",
"% x_ac(3,fmax): Analysis and forecast from AC scheme, with first\n",
"% index representing x, y or z variable\n",
"% x_oi(3,fmax): Analysis and forecast from OI scheme, with first\n",
"% index representing x, y or z variable\n",
"% x_kf(3,fmax): Analysis and forecast from Kalman Filter, with first\n",
"% index representing x, y or z variable\n",
"%\n",
"%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
"\n"
],
"metadata": {
"id": "9HpgUAWU0-Ot"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"## 1. INPUTS"
],
"metadata": {
"id": "mIYDfN2UA5Xh"
}
},
{
"cell_type": "code",
"source": [
"%%octave\n",
"\n",
"%---------------------------------------------------------------------\n",
"%%%%%%%%%%1. INPUTS\n",
"%---------------------------------------------------------------------\n",
"\n",
"%----- Please choose an assimilation scheme -----\n",
"% 1) Successive correction\n",
"% 2) Optimal Interpolation\n",
"% 3) Analysis Correction\n",
"% 4) Kalman Filter\n",
"% Select a menu number:\n",
"l_assim = 1;\n",
"l_assim=l_assim-1;\n",
"\n",
"if (l_assim==0 | l_assim==2)\n",
" %----- How many iterations? -----\n",
" % 1\n",
" % 2\n",
" % 3\n",
" % 4\n",
" % 5\n",
" n_its = 4\n",
"end\n",
"if l_assim~=0\n",
"\n",
"%----- Use correct weighting matrices? -----\n",
"% 1) No\n",
"% 2) Yes\n",
" l_weight = 2;\n",
" l_weight=l_weight-1;\n",
"end\n",
"fmax=1000;\t%time of assimilation + forecast\n",
"q=0:0.01:fmax*0.01-0.01;"
],
"metadata": {
"id": "XZiEmKcTLToq"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"## 2. TRUE SOLUTION\n",
"### True solution of lorenz equations using 2nd order Runge-Kutta method"
],
"metadata": {
"id": "OZ33hAOVBD-F"
}
},
{
"cell_type": "code",
"source": [
"%%octave\n",
"\n",
"%-------------------------------------------------------------------------\n",
"%%%%%%%%%2. TRUE SOLUTION\n",
"%-------------------------------------------------------------------------\n",
"%true solution of lorenz equations using 2nd order Runge-Kutta method\n",
"x=zeros(fmax,1);\n",
"y=zeros(fmax,1);\n",
"z=zeros(fmax,1);"
],
"metadata": {
"id": "bmr1e2-Q7Clz"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"### parameters"
],
"metadata": {
"id": "bGlX8KlFBn1k"
}
},
{
"cell_type": "code",
"source": [
"%%octave\n",
"\n",
"%parameters\n",
"h=0.01; % time step for Runge-Kutta scheme\n",
"s=10.0; % sigma\n",
"r=28; % rho\n",
"b=8/3; % beta"
],
"metadata": {
"id": "eXP_2ayaBTWc"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"### initial contidions"
],
"metadata": {
"id": "CK6XUeR2Bwqk"
}
},
{
"cell_type": "code",
"source": [
"%%octave\n",
"\n",
"%initial conditions\n",
"x(1)=-5.4458d0;\n",
"y(1)=-5.4841d0;\n",
"z(1)=22.5606d0;\n",
"\n",
"for i=1:fmax-1\n",
"k1x(i)=s*y(i)-s*x(i);\n",
"k1y(i)=r*x(i)-y(i)-x(i)*z(i);\n",
"k1z(i)=x(i)*y(i)-b*z(i);\n",
"\n",
"k2x(i)=s*(y(i)+h*k1y(i))-s*(x(i)+h*k1x(i));\n",
"k2y(i)=r*(x(i)+h*k1x(i))-(y(i)+h*k1y(i))-(x(i)+h*k1x(i))*(z(i)+h*k1z(i));\n",
"k2z(i)=(x(i)+h*k1x(i))*(y(i)+h*k1y(i))-b*(z(i)+h*k1z(i));\n",
"\n",
"x(i+1)=x(i)+0.5*h*(k1x(i)+k2x(i));\n",
"y(i+1)=y(i)+0.5*h*(k1y(i)+k2y(i));\n",
"z(i+1)=z(i)+0.5*h*(k1z(i)+k2z(i));\n",
"\n",
"end\n",
"disp(['* done true solution'])"
],
"metadata": {
"id": "eGauO1ooBshO"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"### Plot true solution of x"
],
"metadata": {
"id": "Iniw4SxLB4pl"
}
},
{
"cell_type": "code",
"source": [
"%%octave\n",
"\n",
"%Plot true solution of x\n",
"h1=figure(1);\n",
"clf;\n",
"plot(q,x,'k--');\n",
"legend(\"Truth\")\n",
"xlabel(\"Time\")\n",
"ylabel(\"x\")\n",
"title(\"Solution for x\")\n",
"hold on"
],
"metadata": {
"id": "eiP4BAus7SKN"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"### Plot true solution of z"
],
"metadata": {
"id": "HFX7s8z-CBRG"
}
},
{
"cell_type": "code",
"source": [
"%%octave\n",
"\n",
"h2=figure(2);\n",
"clf;\n",
"plot(q,z,'k--');\n",
"legend(\"Truth\")\n",
"xlabel(\"Time\")\n",
"ylabel(\"z\")\n",
"title('Solution for z')\n",
"hold on"
],
"metadata": {
"id": "VEmxyEFX8eyR"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"## 3. BACKGROUND SOLUTION\n",
"### Produce a background solution which starts from some perturbed initial conditions"
],
"metadata": {
"id": "e1lUNJMbCG87"
}
},
{
"cell_type": "code",
"source": [
"%%octave\n",
"\n",
"%-----------------------------------------------------------------------\n",
"%%%%%%%%%3. BACKGROUND SOLUTION\n",
"%-----------------------------------------------------------------------\n",
"%Produce a background solution which starts from some perturbed initial\n",
"%conditions.\n",
"\n",
"xb=zeros(fmax,1);\n",
"yb=zeros(fmax,1);\n",
"zb=zeros(fmax,1);"
],
"metadata": {
"id": "Q32Bz58r7d8K"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"### Background solution initial condition"
],
"metadata": {
"id": "Ixn-ri7YDJ1j"
}
},
{
"cell_type": "code",
"source": [
"%%octave\n",
"\n",
"%Background solution initial conditions\n",
"%bg_noise=randn(3,1);\n",
"%xb(1)=x(1)+bg_noise(1)*2.0d0;\n",
"%yb(1)=y(1)+bg_noise(2)*2.0d0;\n",
"%zb(1)=z(1)+bg_noise(3)*2.0d0;\n",
"xb(1)=-5.9d0;\n",
"yb(1)=-5.0d0;\n",
"zb(1)=24.0d0;"
],
"metadata": {
"id": "KKDuSi3XDDY7"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"### Loop over time"
],
"metadata": {
"id": "Me3ZlKoRDjiT"
}
},
{
"cell_type": "code",
"source": [
"%%octave\n",
"\n",
"%start loop over time\n",
"for i=1:fmax-1\n",
"\n",
"k1xb(i)=s*yb(i)-s*xb(i);\n",
"k1yb(i)=r*xb(i)-yb(i)-xb(i)*zb(i);\n",
"k1zb(i)=xb(i)*yb(i)-b*zb(i);\n",
"\n",
"k2xb(i)=s*(yb(i)+h*k1yb(i))-s*(xb(i)+h*k1xb(i));\n",
"k2yb(i)=r*(xb(i)+h*k1xb(i))-(yb(i)+h*k1yb(i))-(xb(i)+h*k1xb(i))...\n",
"*(zb(i)+h*k1zb(i));\n",
"k2zb(i)=(xb(i)+h*k1xb(i))*(yb(i)+h*k1yb(i))-b*(zb(i)+h*k1zb(i));\n",
"\n",
"xb(i+1)=xb(i)+0.5*h*(k1xb(i)+k2xb(i));\n",
"yb(i+1)=yb(i)+0.5*h*(k1yb(i)+k2yb(i));\n",
"zb(i+1)=zb(i)+0.5*h*(k1zb(i)+k2zb(i));\n",
"\n",
"end\n",
"%end loop over time\n",
"disp(['* done background solution'])"
],
"metadata": {
"id": "Q4yG11TWDbuS"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"### Plot background solution"
],
"metadata": {
"id": "HKutT7p1ESh8"
}
},
{
"cell_type": "code",
"source": [
"%%octave\n",
"\n",
"%Plot background solution\n",
"figure(h1)\n",
"plot(q,x,'k--');\n",
"hold on\n",
"plot(q,xb,'b-.');\n",
"legend(\"Truth\",\"Background\")\n",
"hold on"
],
"metadata": {
"id": "YgPQ6H7eEEfK"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"source": [
"%%octave\n",
"figure(h2)\n",
"plot(q,z,'k--');\n",
"hold on\n",
"plot(q,zb,'b-.');\n",
"legend(\"Truth\",\"Background\")\n",
"hold on\n",
"\n"
],
"metadata": {
"id": "0_H0bBL9G6pP"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"## 4. OBSERVATIONS"
],
"metadata": {
"id": "X2bPqh7MEh2c"
}
},
{
"cell_type": "markdown",
"source": [
"### Choose number of time steps between observations"
],
"metadata": {
"id": "apQJBpVywEew"
}
},
{
"cell_type": "code",
"source": [
"%%octave\n",
"\n",
"%-----------------------------------------------------------------------------------\n",
"%%%%%%%4 OBSERVATIONS\n",
"%-----------------------------------------------------------------------------------\n",
"\n",
"%----- How many time steps between observations? -----\n",
"% 25\n",
"% 50\n",
"% 100\n",
"% 200\n",
"ob_f=25\n",
"\n",
"%----- Noise on observations? -----\n",
"% 1) No\n",
"% 2) Yes\n",
"% Select a menu number:\n",
"l_noise=2;\n",
"l_noise=l_noise-1;\n",
"tmax=600+ob_f; %time of assimilation (janela de assimilação)"
],
"metadata": {
"id": "zx4ipHYUEcaj"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"### The observations are calculated from the true solution at certain points with noise added."
],
"metadata": {
"id": "-uKbAaNWwlKQ"
}
},
{
"cell_type": "code",
"source": [
"%%octave\n",
"\n",
"%The observations are calculated from the true solution at certain\n",
"%points with noise added.\n",
"\n",
"xob=zeros(tmax,1);\n",
"yob=zeros(tmax,1);\n",
"zob=zeros(tmax,1);\n",
"\n",
"nobs=tmax/ob_f;\t%number of observations\n",
"vec=1:ob_f:tmax;\n",
"\n",
"if l_noise==1\n",
" l_read_noise=1;\n",
" l_read_noise=l_read_noise-1;\n",
" if (l_read_noise==0)\n",
" sc_x_noise=randn(nobs,1);\n",
" sc_y_noise=randn(nobs,1);\n",
" sc_z_noise=randn(nobs,1);\n",
" save ('noise.mat','sc_x_noise','sc_y_noise','sc_z_noise');\n",
" elseif (l_read_noise==1)\n",
" load('noise.mat','sc_x_noise','sc_y_noise','sc_z_noise');\n",
" else\n",
" error('Invalid input')\n",
" end\n",
" sd=0.1;\n",
" var=sqrt(sd);\n",
" RX=var*sc_x_noise;\n",
" RY=var*sc_y_noise;\n",
" RZ=var*sc_z_noise;\n",
"else\n",
" for i=1:nobs;\n",
" RX(i)=0.0;\n",
" RY(i)=0.0;\n",
" RZ(i)=0.0;\n",
" end\n",
"end\n",
"\n",
"j=1;\n",
"for i=vec\n",
"xob(i)=x(i)+RX(j);\n",
"yob(i)=y(i)+RY(j);\n",
"zob(i)=z(i)+RZ(j);\n",
"j=j+1;\n",
"end\n",
"disp(['* done observations'])\n",
"\n"
],
"metadata": {
"id": "8KhPtHnXwV-q"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"### Plot x observations"
],
"metadata": {
"id": "m1-TAPgJ7INl"
}
},
{
"cell_type": "code",
"source": [
"%%octave\n",
"\n",
"%Plot x observations\n",
"v=0.01*vec;\n",
"figure(h1)\n",
"plot(q,x,'k--');\n",
"hold on\n",
"plot(q,xb,'b-.');\n",
"plot(v(2:nobs),xob(vec(2:nobs)),'om');\n",
"legend(\"Truth\",\"Background\",\"Observations\")\n",
"hold on\n",
"\n"
],
"metadata": {
"id": "LMuSMG9533bF"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"### Plot z observations"
],
"metadata": {
"id": "YHKwdGUn7qb-"
}
},
{
"cell_type": "code",
"source": [
"%%octave\n",
"\n",
"figure(h2)\n",
"plot(q,z,'k--');\n",
"hold on\n",
"plot(q,zb,'b-.');\n",
"plot(v(2:nobs),zob(vec(2:nobs)),'om');\n",
"legend(\"Truth\",\"Background\",\"Observations\")\n",
"hold on\n",
"\n"
],
"metadata": {
"id": "hb4lbcBm7P10"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"## 5. ERROR COVARIANCE MATRICES"
],
"metadata": {
"id": "8uGTh48s7aVz"
}
},
{
"cell_type": "code",
"source": [
"%%octave\n",
"\n",
"%-----------------------------------------------------------------------------------\n",
"%%%%%%%5. ERROR COVARIANCE MATRICES\n",
"%-----------------------------------------------------------------------------------\n",
"% l_test=menu_asl('Click to perform analysis','OK');\n",
"if l_assim~=0\n",
"\n",
"BX=zeros(3,3);\n",
"BXi=zeros(3,3,tmax);\n",
"R=zeros(3,3);\n",
"Ri=zeros(3,3,nobs);\n",
"\n",
"%Observation error covariance matrix, R.\n",
"%Calculated to be the average errors over time between\n",
"%the observations and true solution at the observation points.\n",
"R=zeros(3,3);\n",
" for i=vec\n",
" Ri(1,1,i)=(x(i)-xob(i))'*(x(i)-xob(i));\n",
" Ri(1,2,i)=(x(i)-xob(i))'*(y(i)-yob(i));\n",
" Ri(1,3,i)=(x(i)-xob(i))'*(z(i)-zob(i));\n",
" Ri(2,1,i)=(y(i)-yob(i))'*(x(i)-xob(i));\n",
" Ri(2,2,i)=(y(i)-yob(i))'*(y(i)-yob(i));\n",
" Ri(2,3,i)=(y(i)-yob(i))'*(z(i)-zob(i));\n",
" Ri(3,1,i)=(z(i)-zob(i))'*(x(i)-xob(i));\n",
" Ri(3,2,i)=(z(i)-zob(i))'*(y(i)-yob(i));\n",
" Ri(3,3,i)=(z(i)-zob(i))'*(z(i)-zob(i));\n",
" end\n",
"\n",
"if l_weight==1\n",
" Rj=zeros(3,3,nobs);\n",
" for j=1:3\n",
" for k=1:3\n",
" Rj(j,k,:)=Ri(j,k,vec);\n",
" R(j,k)=(1/nobs)*sum(Rj(j,k,:));\n",
" end\n",
" end\n",
"else\n",
" R=eye(3,3); % Set to identity matrix for l_weight==0\n",
"end\n",
"\n",
"%Background error covariance matrix.\n",
"%Calculated to be the average errors over time between the\n",
"%background and true solutions.\n",
"B=zeros(3,3);\n",
"if l_weight==1\n",
" for i=1:tmax\n",
" BXi(1,1,i)=(x(i)-xb(i))'*(x(i)-xb(i));\n",
" BXi(1,2,i)=(x(i)-xb(i))'*(y(i)-yb(i));\n",
" BXi(1,3,i)=(x(i)-xb(i))'*(z(i)-zb(i));\n",
" BXi(2,1,i)=(y(i)-yb(i))'*(x(i)-xb(i));\n",
" BXi(2,2,i)=(y(i)-yb(i))'*(y(i)-yb(i));\n",
" BXi(2,3,i)=(y(i)-yb(i))'*(z(i)-zb(i));\n",
" BXi(3,1,i)=(z(i)-zb(i))'*(x(i)-xb(i));\n",
" BXi(3,2,i)=(z(i)-zb(i))'*(y(i)-yb(i));\n",
" BXi(3,3,i)=(z(i)-zb(i))'*(z(i)-zb(i));\n",
" end\n",
" for j=1:3\n",
" for k=1:3\n",
" BX(j,k)=(1/tmax)*sum(BXi(j,k,:));\n",
" end\n",
" end\n",
"else\n",
" BX=eye(3,3); % Set to identity matrix for l_weight==0\n",
"end % l_weight==1\n",
"\n",
"end % l_assim~=0"
],
"metadata": {
"id": "3KcRxVr97jok"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"### 6. DATA ASSIMILATION"
],
"metadata": {
"id": "iEyiPjj7eK9u"
}
},
{
"cell_type": "code",
"source": [
"%%octave\n",
"\n",
"%-----------------------------------------------------------------------------------\n",
"%%%%%%%%%%6. DATA ASSIMILATION\n",
"%-----------------------------------------------------------------------------------\n",
"\n",
"\n",
"if l_assim==0\n",
"%%%%%%%%%%%%%%%%%% SUCCESSIVE CORRECTIONS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
"disp('Performing successive correction analysis')\n",
"assim_str='Successive correction';\n",
"x_ob=zeros(3,fmax);\n",
"x_sc=zeros(3,fmax);\n",
"WX=zeros(3,3);\n",
"\n",
"%state weighting matrix\n",
"WX=0.5*eye(3,3);\n",
"\n",
"%forecast at the initial time is the background there.\n",
"xfc(1)=xb(1);\n",
"yfc(1)=yb(1);\n",
"zfc(1)=zb(1);\n",
"\n",
"x_ob=[xob';yob';zob'];\n",
"x_sc(:,1)=[xb(1);yb(1);zb(1)];\n",
"\n",
"\n",
"%Start Loop over time\n",
"for i=1:fmax-1\n",
"\n",
" if i>=tmax\n",
" x_ob(:,i+1)=0;\n",
" end\n",
"\n",
"%Forecast to time step i+1\n",
" xfc(i)=x_sc(1,i);\n",
" yfc(i)=x_sc(2,i);\n",
" zfc(i)=x_sc(3,i);\n",
"\n",
" k1xfc(i)=s*yfc(i)-s*xfc(i);\n",
" k1yfc(i)=r*xfc(i)-yfc(i)-xfc(i)*zfc(i);\n",
" k1zfc(i)=xfc(i)*yfc(i)-b*zfc(i);\n",
"\n",
" k2xfc(i)=s*(yfc(i)+h*k1yfc(i))-s*(xfc(i)+h*k1xfc(i));\n",
" k2yfc(i)=r*(xfc(i)+h*k1xfc(i))-(yfc(i)+h*k1yfc(i))-(xfc(i)+h*k1xfc(i))...\n",
" *(zfc(i)+h*k1zfc(i));\n",
" k2zfc(i)=(xfc(i)+h*k1xfc(i))*(yfc(i)+h*k1yfc(i))-b*(zfc(i)+h*k1zfc(i));\n",
"\n",
" xfc(i+1)=xfc(i)+0.5*h*(k1xfc(i)+k2xfc(i));\n",
" yfc(i+1)=yfc(i)+0.5*h*(k1yfc(i)+k2yfc(i));\n",
" zfc(i+1)=zfc(i)+0.5*h*(k1zfc(i)+k2zfc(i));\n",
"\n",
"%If there is an observation then produce an analysis\n",
" if x_ob(:,i+1)~=0\n",
" x_sc(:,i+1)=[xfc(i+1);yfc(i+1);zfc(i+1)];\n",
"\n",
" for j=1:n_its\n",
" x_sc(:,i+1)=x_sc(:,i+1)+WX*(x_ob(:,i+1)-x_sc(:,i+1));\n",
" end\n",
"\n",
" else\n",
"\n",
" x_sc(:,i+1)=[xfc(i+1);yfc(i+1);zfc(i+1)];\n",
"\n",
" end\n",
"end\n",
"\n",
"\n",
"\n",
"elseif l_assim==1\n",
"%%%%%%%%%%%%%%%%%% OPTIMAL INTERPOLATION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
"disp('Performing OI analysis')\n",
"assim_str='Optimal interpolation';\n",
"\n",
"x_ob=zeros(3,fmax);\n",
"x_oi=zeros(3,fmax);\n",
"WX=zeros(3,3);\n",
"\n",
"xfc(1)=xb(1);\n",
"yfc(1)=yb(1);\n",
"zfc(1)=zb(1);\n",
"\n",
"%state weighting matrix\n",
"WX=BX*pinv(BX+R);\n",
"\n",
"x_ob=[xob';yob';zob'];\n",
"x_oi(:,1)=[xb(1);yb(1);zb(1)];\n",
"\n",
"%Start loop over time\n",
"for i=1:fmax-1\n",
"\n",
" if i>=tmax\n",
" x_ob(:,i+1)=0;\n",
" end\n",
"\n",
"%Forecast to time step i+1\n",
" xfc(i)=x_oi(1,i);\n",
" yfc(i)=x_oi(2,i);\n",
" zfc(i)=x_oi(3,i);\n",
"\n",
" k1xfc(i)=s*yfc(i)-s*xfc(i);\n",
" k1yfc(i)=r*xfc(i)-yfc(i)-xfc(i)*zfc(i);\n",
" k1zfc(i)=xfc(i)*yfc(i)-b*zfc(i);\n",
"\n",
" k2xfc(i)=s*(yfc(i)+h*k1yfc(i))-s*(xfc(i)+h*k1xfc(i));\n",
" k2yfc(i)=r*(xfc(i)+h*k1xfc(i))-(yfc(i)+h*k1yfc(i))-(xfc(i)+h*k1xfc(i))...\n",
" *(zfc(i)+h*k1zfc(i));\n",
" k2zfc(i)=(xfc(i)+h*k1xfc(i))*(yfc(i)+h*k1yfc(i))-b*(zfc(i)+h*k1zfc(i));\n",
"\n",
" xfc(i+1)=xfc(i)+0.5*h*(k1xfc(i)+k2xfc(i));\n",
" yfc(i+1)=yfc(i)+0.5*h*(k1yfc(i)+k2yfc(i));\n",
" zfc(i+1)=zfc(i)+0.5*h*(k1zfc(i)+k2zfc(i));\n",
"\n",
"%If there is an observation then produce analysis\n",
" if x_ob(:,i+1)~=0\n",
"\n",
" x_oi(:,i+1)=[xfc(i+1);yfc(i+1);zfc(i+1)];\n",
" x_oi(:,i+1)=x_oi(:,i+1)+WX*(x_ob(:,i+1)-x_oi(:,i+1));\n",
"\n",
" else\n",
"\n",
" x_oi(:,i+1)=[xfc(i+1);yfc(i+1);zfc(i+1)];\n",
"\n",
" end\n",
"end\n",
"\n",
"elseif l_assim==2\n",
"%%%%%%%%%%%%%%%%%%%%%%%%%% ANALYSIS CORRECTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
"disp('Performing AC analysis')\n",
"assim_str='Analysis correction';\n",
"\n",
"if l_noise==0\n",
" R=eye(3,3);\n",
"end\n",
"x_ob=zeros(3,fmax);\n",
"x_ac=zeros(3,fmax);\n",
"y_ac=zeros(3,fmax);\n",
"\n",
"WX=zeros(3,3); %=BXR^{-1}\n",
"QX=zeros(3,3); %=(WX+I)^{-1}\n",
"\n",
"xfc(1)=xb(1);\n",
"yfc(1)=yb(1);\n",
"zfc(1)=zb(1);\n",
"\n",
"%state weighting matrix\n",
"WX=BX*pinv(R)\n",
"QX=pinv(WX+eye(3,3))\n",
"\n",
"x_ob=[xob';yob';zob'];\n",
"x_ac(:,1)=[xb(1);yb(1);zb(1)];\n",
"y_ac(:,1)=x_ob(:,1);\n",
"\n",
"%Start loop over time\n",
"for i=1:fmax-1\n",
"\n",
" if i>=tmax\n",
" x_ob(:,i+1)=0;\n",
" end\n",
"\n",
"%Forecast to time step i+1\n",
" xfc(i)=x_ac(1,i);\n",
" yfc(i)=x_ac(2,i);\n",
" zfc(i)=x_ac(3,i);\n",
"\n",
" k1xfc(i)=s*yfc(i)-s*xfc(i);\n",
" k1yfc(i)=r*xfc(i)-yfc(i)-xfc(i)*zfc(i);\n",
" k1zfc(i)=xfc(i)*yfc(i)-b*zfc(i);\n",
"\n",
" k2xfc(i)=s*(yfc(i)+h*k1yfc(i))-s*(xfc(i)+h*k1xfc(i));\n",
" k2yfc(i)=r*(xfc(i)+h*k1xfc(i))-(yfc(i)+h*k1yfc(i))-(xfc(i)+h*k1xfc(i))...\n",
" *(zfc(i)+h*k1zfc(i));\n",
" k2zfc(i)=(xfc(i)+h*k1xfc(i))*(yfc(i)+h*k1yfc(i))-b*(zfc(i)+h*k1zfc(i));\n",
"\n",
" xfc(i+1)=xfc(i)+0.5*h*(k1xfc(i)+k2xfc(i));\n",
" yfc(i+1)=yfc(i)+0.5*h*(k1yfc(i)+k2yfc(i));\n",
" zfc(i+1)=zfc(i)+0.5*h*(k1zfc(i)+k2zfc(i));\n",
"\n",
"%If there is an observation then produce an analysis\n",
" if x_ob(:,i+1)~=0\n",
" x_ac(:,i+1)=[xfc(i+1);yfc(i+1);zfc(i+1)];\n",
" y_ac(:,i+1)=x_ob(:,i+1);\n",
"\n",
" for j=1:n_its\n",
" x_ac(:,i+1)=x_ac(:,i+1)+WX*QX*(y_ac(:,i+1)-x_ac(:,i+1));\n",
" y_ac(:,i+1)=y_ac(:,i+1)-QX*(y_ac(:,i+1)-x_ac(:,i+1));\n",
" end\n",
"\n",
" else\n",
"\n",
" x_ac(:,i+1)=[xfc(i+1);yfc(i+1);zfc(i+1)];\n",
"\n",
" end\n",
"end\n",
"%End loop over time\n",
"\n",
"elseif l_assim==3\n",
"%%%%%%%%%%%%%%%%%%%%%%%%%% KALMAN FILTER %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n",
"disp('Performing kalman filter analysis')\n",
"assim_str='Kalman filter';\n",
"\n",
"KX=zeros(3,3,fmax);\n",
"Pfx=zeros(3,3,fmax);\n",
"Pax=zeros(3,3,fmax);\n",
"xfc=zeros(tmax);\n",
"yfc=zeros(tmax);\n",
"zfc=zeros(tmax);\n",
"x_kf=zeros(3,fmax);\n",
"x_fc=zeros(3,fmax);\n",
"x_ob=zeros(3,fmax);\n",
"Mx=zeros(3,3,fmax);\n",
"Qx=zeros(3,3,fmax);\n",
"\n",
"%put obs into one array\n",
"x_ob=[xob';yob';zob'];\n",
"\n",
"%set forecast at i=1 to be the background guess there\n",
"xfc(1)=xb(1);\n",
"yfc(1)=yb(1);\n",
"zfc(1)=zb(1);\n",
"\n",
"%put forecast values into one array\n",
"x_fc(:,1)=[xb(1);yb(1);zb(1)];\n",
"\n",
"\n",
"if l_weight==1\n",
"%forecast error covariance at i=1 is the background error\n",
"%covariance there\n",
" Pfx(:,:,1)=BXi(:,:,1);\n",
"else\n",
"%begin with some incorrect guess\n",
" Pfx(:,:,1)=eye(3);\n",
"end\n",
"\n",
"%calculate the gain matrix at i=1\n",
"KX(:,:,1)=Pfx(:,:,1)*pinv(Pfx(:,:,1)+Ri(:,:,1));\n",
"\n",
"x_kf(:,1)=x_fc(:,1);\n",
"\n",
"%calculate analysis error covariance matrix at i=1\n",
"Pax(:,:,1)=(eye(3,3)-KX(:,:,1))*Pfx(:,:,1);\n",
"\n",
"\n",
"for i=1:fmax-1\n",
" if i>=tmax\n",
" x_ob(:,i+1)=0;\n",
" end\n",
"\n",
" xfc(i)=x_kf(1,i);\n",
" yfc(i)=x_kf(2,i);\n",
" zfc(i)=x_kf(3,i);\n",
"\n",
"%tangent-linear model - linearised about the forecast state\n",
" Mx(1,1,i)=1.0-(h*s)+(((h*h*s)/2.0)*(r+s-zfc(i)));\n",
" Mx(1,2,i)=h*s-(((h*h*s)/2.0)*(1.0+s));\n",
" Mx(1,3,i)=-(h*h*s*xfc(i))/2.0;\n",
" Mx(2,1,i)=(h*(r-zfc(i)))+(((h*h)/2.0)*(zfc(i)-(r*s)-r+(b*zfc(i))...\n",
" +(s*zfc(i))-(2.0*xfc(i)*yfc(i))))+(((h*h*h*s)/2.0)*((b*zfc(i))-...\n",
" (yfc(i)*yfc(i))+(2.0*xfc(i)*yfc(i))));\n",
" Mx(2,2,i)=1.0-h+(((h*h)/2.0)*((s*r)+1.0-(xfc(i)*xfc(i))-(s*zfc(i))))...\n",
" +(((h*h*h*s)/2.0)*((b*zfc(i))-(2.0*xfc(i)*yfc(i))+(xfc(i)*xfc(i))));\n",
" Mx(2,3,i)=(-h*xfc(i))+(((h*h)/2.0)*(xfc(i)+(b*xfc(i))-(s*yfc(i))+...\n",
" (s*xfc(i))))+(((h*h*h*s*b)/2.0)*(yfc(i)+xfc(i)));\n",
" Mx(3,1,i)=(h*yfc(i))+(((h*h)/2.0)*((2.0*r*xfc(i))-yfc(i)-(2.0*xfc(i)...\n",
" *zfc(i))-(s*yfc(i))-(b*yfc(i))))+(((h*h*h*s)/2.0)*((r*yfc(i))-(yfc(i)...\n",
" *zfc(i))-(2.0*r*xfc(i))+yfc(i)+(2.0*xfc(i)*zfc(i))));\n",
" Mx(3,2,i)=(h*xfc(i))+(((h*h)/2.0)*((s*yfc(i))-xfc(i)-(s*xfc(i))-...\n",
" (b*xfc(i))))+(((h*h*h*s)/2.0)*((r*xfc(i))-(2.0*yfc(i))-(xfc(i)*zfc(i))...\n",
" +xfc(i)));\n",
" Mx(3,3,i)=1.0-(h*b)+(((h*h)/2.0)*((b*b)-(xfc(i)*xfc(i))))+...\n",
" (((h*h*h*s)/2.0)*((xfc(i)*xfc(i))-(xfc(i)*yfc(i))));\n",
"\n",
"%Choose something for the model error covariance matrix.\n",
"%These are from Evensen, 1997.\n",
" Qx(1,1,i)=0.1491;\n",
" Qx(1,2,i)=0.1505;\n",
" Qx(1,3,i)=0.0007;\n",
" Qx(2,1,i)=0.1505;\n",
" Qx(2,2,i)=0.9048;\n",
" Qx(2,3,i)=0.0014;\n",
" Qx(3,1,i)=0.0007;\n",
" Qx(3,2,i)=0.0014;\n",
" Qx(3,3,i)=0.9180;\n",
"\n",
"%forecast the state\n",
" k1xfc(i)=s*yfc(i)-s*xfc(i);\n",
" k1yfc(i)=r*xfc(i)-yfc(i)-xfc(i)*zfc(i);\n",
" k1zfc(i)=xfc(i)*yfc(i)-b*zfc(i);\n",
"\n",
" k2xfc(i)=s*(yfc(i)+h*k1yfc(i))-s*(xfc(i)+h*k1xfc(i));\n",
" k2yfc(i)=r*(xfc(i)+h*k1xfc(i))-(yfc(i)+h*k1yfc(i))-(xfc(i)+h*k1xfc(i))...\n",
" *(zfc(i)+h*k1zfc(i));\n",
" k2zfc(i)=(xfc(i)+h*k1xfc(i))*(yfc(i)+h*k1yfc(i))-b*(zfc(i)+h*k1zfc(i));\n",
"\n",
" xfc(i+1)=xfc(i)+0.5*h*(k1xfc(i)+k2xfc(i));\n",
" yfc(i+1)=yfc(i)+0.5*h*(k1yfc(i)+k2yfc(i));\n",
" zfc(i+1)=zfc(i)+0.5*h*(k1zfc(i)+k2zfc(i));\n",
"\n",
" x_fc(:,i+1)=[xfc(i+1);yfc(i+1);zfc(i+1)];\n",
"\n",
"%forecast error covariance matrix\n",
" Pfx(:,:,i+1)=Mx(:,:,i)*Pax(:,:,i)*Mx(:,:,i)'+Qx(:,:,i);\n",
"\n",
"%If there is an observation at this time\n",
" if x_ob(:,i+1)~=0\n",
"\n",
"%calculate the gain matrices\n",
" KX(:,:,i+1)=Pfx(:,:,i+1)*pinv(Pfx(:,:,i+1)+Ri(:,:,i+1));\n",
"\n",
"%produce an analysis\n",
" x_kf(:,i+1)=x_fc(:,i+1)+KX(:,:,i+1)*(x_ob(:,i+1)-x_fc(:,i+1));\n",
"\n",
"%calculate analysis error covariance matrices\n",
" Pax(:,:,i+1)=(eye(3,3)-KX(:,:,i+1))*Pfx(:,:,i+1);\n",
"\n",
" else\n",
"\n",
" x_kf(:,i+1)=x_fc(:,i+1);\n",
" Pax(:,:,i+1)=Pfx(:,:,i+1);\n",
"\n",
" end\n",
"\n",
"end\n",
"%End loop over time\n",
"\n",
"end\n",
"%End if on type of assimilation\n",
"\n",
"disp(['*** FINISHED ***'])"
],
"metadata": {
"id": "UNUdvYKneCbe"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"## 7. PLOT RESULTS"
],
"metadata": {
"id": "nyNs4Q-O2Qry"
}
},
{
"cell_type": "markdown",
"source": [
"### Plot analysed solution of x"
],
"metadata": {
"id": "yASAkISp2a89"
}
},
{
"cell_type": "code",
"source": [
"%%octave\n",
"\n",
"%-----------------------------------------------------------------------\n",
"%%%%%%%%%%7. PLOT RESULTS\n",
"%-----------------------------------------------------------------------\n",
"\n",
"%Plot analysed solution of x\n",
"figure(h1)\n",
"plot(q,x,'k--');\n",
"hold on\n",
"plot(q,xb,'b-.');\n",
"plot(v(2:nobs),xob(vec(2:nobs)),'om');\n",
"if l_assim==0\n",
" plot(q,x_sc(1,:),'r-');\n",
"elseif l_assim==1\n",
" plot(q,x_oi(1,:),'r-');\n",
"elseif l_assim==2\n",
" plot(q,x_ac(1,:),'r-');\n",
"elseif l_assim==3\n",
" plot(q,x_kf(1,:),'r-');\n",
"end\n",
"legend(\"Truth\",\"Background\",\"Observations\",\"Analysis\")\n",
"yminmax=get(gca,'YLim');\n",
"yspace=(yminmax(2)-yminmax(1))*0.01;\n",
"yvals=yminmax(1):yspace:yminmax(2);\n",
"xval=(tmax-ob_f)*h;\n",
"%%line(xval,yvals,'LineStyle',':','Color','k')\n",
"hold on\n"
],
"metadata": {
"id": "5baYdWyl2M-K"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"### Plot analysed solution of Z"
],
"metadata": {
"id": "ltN_czPf2nXJ"
}
},
{
"cell_type": "code",
"source": [
"%%octave\n",
"\n",
"figure(h2);\n",
"plot(q,z,'k--');\n",
"hold on\n",
"plot(q,zb,'b-.');\n",
"plot(v(2:nobs),zob(vec(2:nobs)),'om');\n",
"%Plot analysed solution of z\n",
"if l_assim==0\n",
" plot(q,x_sc(3,:),'r-');\n",
"elseif l_assim==1\n",
" plot(q,x_oi(3,:),'r-');\n",
"elseif l_assim==2\n",
" plot(q,x_ac(3,:),'r-');\n",
"elseif l_assim==3\n",
" plot(q,x_kf(3,:),'r-');\n",
"end\n",
"legend(\"Truth\",\"Background\",\"Observations\",\"Analysis\")\n",
"yminmax=get(gca,'YLim');\n",
"yspace=(yminmax(2)-yminmax(1))*0.01;\n",
"yvals=yminmax(1):yspace:yminmax(2);\n",
"xval=(tmax-ob_f)*h;\n",
"%line(xval,yvals,'LineStyle',':','Color','k')\n",
"hold on"
],
"metadata": {
"id": "65TQqQLl2kV6"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [],
"metadata": {
"id": "tS_b6MDs29RO"
}
},
{
"cell_type": "code",
"source": [
"%%octave\n",
"\n",
"if l_assim==0\n",
" x_err = x'-x_sc(1,:);\n",
" z_err = z'-x_sc(3,:);\n",
"elseif l_assim==1\n",
" x_err = x'-x_oi(1,:);\n",
" z_err = z'-x_oi(3,:);\n",
"elseif l_assim==2\n",
" x_err = x'-x_ac(1,:);\n",
" z_err = z'-x_ac(3,:);\n",
"elseif l_assim==3\n",
" x_err = x'-x_kf(1,:);\n",
" z_err = z'-x_kf(3,:);\n",
"end\n",
"\n",
"h3=figure(3);\n",
"clf;\n",
"subplot(2,1,1)\n",
"plot(q,x_err,'k-');\n",
"xlabel('Time');ylabel('Error')\n",
"title('Plot of error (Truth-Analysis) in x against time')\n",
"yminmax=get(gca,'YLim');\n",
"yspace=(yminmax(2)-yminmax(1))*0.01;\n",
"yvals=yminmax(1):yspace:yminmax(2);\n",
"%%line(xval,yvals,'LineStyle',':','Color','k')\n",
"%\n",
"subplot(2,1,2)\n",
"plot(q,z_err,'k-')\n",
"xlabel('Time');ylabel('Error')\n",
"title('Plot of error (Truth-Analysis) in z against time')\n",
"yminmax=get(gca,'YLim');\n",
"yspace=(yminmax(2)-yminmax(1))*0.01;\n",
"yvals=yminmax(1):yspace:yminmax(2);\n",
"%%line(xval,yvals,'LineStyle',':','Color','k')\n",
"\n"
],
"metadata": {
"id": "9HksltdQ27Je"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "markdown",
"source": [
"### Write out options chosen"
],
"metadata": {
"id": "txQqnLrk3gOx"
}
},
{
"cell_type": "code",
"source": [
"%%octave\n",
"\n",
"% Determine screen size and set up figure positions\n",
"ss=get(0,'ScreenSize');\n",
"fig_width=0.55*ss(3);\n",
"fig_height=0.4*ss(4);\n",
"pos_1=[0.02*ss(3),0.55*ss(4),fig_width,fig_height];\n",
"pos_2=[0.02*ss(3),0.05*ss(4),fig_width,fig_height];\n",
"pos_3=[0.6*ss(3),0.45*ss(4),0.4*ss(3),0.45*ss(4)];\n",
"pl=0.6*ss(3);\n",
"pb=0.05*ss(4);\n",
"pw=0.3*ss(3);\n",
"ph=0.3*ss(4);\n",
"pos_4=[pl,pb,pw,ph];\n",
"pos_5=[0.02*pw,0.8*ph,0.95*pw,0.2*ph];\n",
"pos_6=[0.02*pw,0.02*ph,0.95*pw,0.75*ph];"
],
"metadata": {
"id": "gy_T7rtwvqpp"
},
"execution_count": null,
"outputs": []
},
{
"cell_type": "code",
"source": [
"%%octave\n",
"%\n",
"%-----------------------------------------------------------------------\n",
"% Write out options chosen\n",
"%\n",
"text1={'List of options chosen'};\n",
"text2={['Analysis scheme: ' assim_str]};\n",
"if (l_assim==0 | l_assim==2)\n",
" text3={['Number of iterations: ' num2str(n_its)]};\n",
"else\n",
" text3='';\n",
"end\n",
"if l_assim~0\n",
" if l_weight==0\n",
" text4={['Using incorrect weighting matrices']};\n",
" else\n",
" text4={['Using correct weighting matrices']};\n",
" end\n",
"else\n",
" text4='';\n",
"end\n",
"text5={['Time steps between observations: ' num2str(ob_f)]};\n",
"if l_noise==0\n",
" text6={['No noise on observations']};\n",
" text7='';\n",
"else\n",
" text6={['Observations have random noise with variance ' num2str(sd)]};\n",
" if l_read_noise==0\n",
" text7={['Noise generated in program and saved to file']};\n",
" else\n",
" text7={['Noise read in from previously generated file']};\n",
" end\n",
"end\n",
"%\n",
"str1=[text2;text3;text4;text5;text6;text7];\n",
"disp(str1)\n",
"\n",
"%%%%%%%%%% END OF PROGRAM %%%%%%%%%%\n"
],
"metadata": {
"id": "vUWy5LRa3arh"
},
"execution_count": null,
"outputs": []
}
]
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment