{
"cells": [
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"pkg load control\n",
"%plot gnuplot\n",
"%plot --format svg\n"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"function dx = f(x,t)\n",
" dx(1) = x(2) - x(1);\n",
" dx(2) = x(1)*(4-x(2));\n",
"end\n",
"xs = lsode(@f,[1,2],0:0.01:10);"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"function phasors(Z)\n",
" figure\n",
" hold on\n",
" \n",
" % Calculate a suitable axis limit based on phasors\n",
" m = 1.2 * max([abs(real(Z)) abs(imag(Z))]);\n",
" \n",
" % Create a new figure window\n",
" #figure\n",
" #hold on\n",
" \n",
" % Plot phasors one at a time\n",
" quiver(0*real(Z),0*imag(Z),real(Z),imag(Z));\n",
" \n",
" % Draw horizontal and vertical axis lines\n",
" plot([-m m],[0 0])\n",
" plot([0 0],[-m m])\n",
" \n",
" % Set axis limits and aspect ratio\n",
" axis([-m,m,-m,m], 'square')\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Current magnitude (A):\n",
" 0.0029972\n",
"Current phase angle (rad):\n",
" 1.2664\n"
]
},
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%plot gnuplot\n",
"%plot --format svg\n",
"\n",
"% RCphasors.m - An RC circuit example\n",
" \n",
"% create a variable j equal to sqrt(-1)\n",
"j = 0 + 1j;\n",
" \n",
"% frequency in Hz and angular frequency in rad/s\n",
"f = 50;\n",
"w = 2*pi*f;\n",
" \n",
"% supply voltage has magnitude 10V and zero phase angle\n",
"% i.e. at its maximum value at t=0\n",
"Vs = 10;\n",
" \n",
"% Component values\n",
"R = 1e3;\n",
"C = 1e-6;\n",
" \n",
"% impedance of capacitor\n",
"Zc = 1 / (j*w*C);\n",
" \n",
"% total equivalent impedance of series R and C\n",
"Zeq = R + Zc;\n",
" \n",
"% Calculate current\n",
"I = Vs / Zeq;\n",
" \n",
"% Calculate capacitor and resistor voltages\n",
"Vc = I * Zc;\n",
"Vr = I * R;\n",
" \n",
"% Display current magnitude and phase\n",
"disp 'Current magnitude (A):'\n",
"disp (abs(I))\n",
"disp 'Current phase angle (rad):'\n",
"disp (angle(I))\n",
" \n",
"% Plot a phasor diagram (scale I to make it visible)\n",
"phasors([Vs 1000*I]);"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"clear;"
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fs=8000;\n",
"a = [ 1 -0.455938 ];\n",
"b = [ 0.544062 ];\n",
"[H,f] = freqz(b, a, 64, fs);\n",
"cutoff = -3 * ones(64);\n",
"octaveup = -9 * ones(64);\n",
"figure();\n",
"plot(f, 20*log10(abs(H)), f, cutoff, f, octaveup);\n",
"xlabel(\"Frequency (Hz)\");\n",
"ylabel(\"Magnitude (dB)\");"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"fs=8000;\n",
"a = [ 1 -0.414213 ];\n",
"b = [ -0.585786 -0.585786] / 2;\n",
"[H,f] = freqz(b, a, 64, fs);\n",
"cutoff = -3 * ones(64);\n",
"octaveup = -9 * ones(64);\n",
"figure(1);\n",
"plot(f, 20*log10(abs(H)), f, cutoff, f, octaveup);\n",
"xlabel(\"Frequency (Hz)\");\n",
"ylabel(\"Magnitude (dB)\");"
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"[x1, x2] = meshgrid(-.5:0.05:0.5, -.5:.05:.5);\n",
"x1dot = -x1 - 2 *x2 .*x1.^2+x2; %Note the use of .* and .^\n",
"x2dot = -x1-x2;\n",
"quiver(x1,x2,x1dot, x2dot)"
]
},
{
"cell_type": "code",
"execution_count": 57,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"%plot gnuplot\n",
"%plot --format svg\n",
"pkg load control\n",
"f = @(t,Y) [Y(2); -sin(Y(1))];\n",
"y1 = linspace(-2,8,20);\n",
"y2 = linspace(-2,2,20);\n",
"[x,y] = meshgrid(y1,y2);\n",
"size(x);\n",
"size(y);\n",
"u = zeros(size(x));\n",
"v = zeros(size(x));\n",
"t=0;\n",
"for i = 1:numel(x)\n",
" Yprime = f(t,[x(i); y(i)]);\n",
" u(i) = Yprime(1);\n",
" v(i) = Yprime(2);\n",
"end\n",
"figure();\n",
"quiver(x,y,u,v,'r'); \n",
"xlabel('y_1');\n",
"ylabel('y_2');\n",
"axis tight equal;\n",
"#hold on\n",
"#for y20 = [0 0.5 1 1.5 2 2.5]\n",
"# [ts,ys] = lsode(\"f\",0,0.5,2.5);\n",
" # plot(ys(:,1),ys(:,2))\n",
" # plot(ys(1,1),ys(1,2),'bo') % starting point\n",
" # plot(ys(end,1),ys(end,2),'ks') % ending point\n",
"#end\n",
"#hold off\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 52,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Define the right-hand side of the equation:\n",
"function ret=f(x,t); ret=x^2; end;\n",
"# t will be on the interval [0,1]; x(0)=0.5\n",
"# t will be the set of moments of time:\n",
"t=(0:0.1:1)';\n",
"# x will be the values of the function at these moments of time.\n",
"x=lsode('f',0.5,t);\n",
"plot(t,x)"
]
},
{
"cell_type": "code",
"execution_count": 62,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"[x1, x2] = meshgrid(-.5:.05:.5, -.5:.05:.5);\n",
"x1dot = 2 * x2 .* (1 + x1.^2 + x2.^2);\n",
"x2dot = -4 * x1 .* (1 + x1.^2 + x2.^2);\n",
"\n",
"% plot vector field\n",
"quiver(x1, x2, x1dot, x2dot)\n",
"grid\n",
"xlabel('x1');\n",
"ylabel('x2');\n",
"\n",
"% some output\n",
"title(\"Phasen-Portrait eines einfachen Systems\")\n",
"print(\"phaseportrait.pdf\", \"-color\", \"-dpdf\")\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 67,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"timer = 0.24000\r\n"
]
},
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"Tend = 30; u0 = [0;0];\n",
"function curr = Diode(u)\n",
" Rd = 10; us = 0.7;\n",
" if (u>=-us) curr=0;\n",
" else curr=Rd*(u+us);\n",
" endif\n",
"endfunction\n",
"function y = circuit(u,t)\n",
" C1 = 1; C2 = 1;\n",
" y = [-1/C1*(Diode(u(1)-10*sin(t))-Diode(u(2)-u(1)));\n",
" 10*cos(t)-1/C2*Diode(u(2)-u(1))];\n",
"endfunction\n",
"t = linspace(0,Tend,100);\n",
"lsode_options(\"absolute tolerance\",1e-5);\n",
"lsode_options(\"relative tolerance\",1e-5);\n",
"t0 = cputime();\n",
"u = lsode('circuit',u0,t);\n",
"timer = cputime()-t0\n",
"plot(t,u(:,2))\n",
"grid on; xlabel('time'); ylabel('tension');"
]
},
{
"cell_type": "code",
"execution_count": 71,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"error: 'rk45' undefined near line 1 column 8\n",
"timer = 0.012000\n",
"error: 'ode_Runge' undefined near line 1 column 14\n",
"timer = 0.0080000\n",
"error: 'uFix' undefined near line 1 column 11\n",
"error: evaluating argument list element number 2\n",
"warning: legend: plot data is empty; setting key labels has no effect\n",
"warning: called from\n",
" legend at line 378 column 9\n",
"error: legend: subscript indices must be either positive integers less than 2^31 or logicals\n",
"error: called from\n",
" legend at line 420 column 15\n"
]
},
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"Tend = 30; u0 = [0;0];\n",
"function curr = Diode(u)\n",
" Rd = 10; us = 0.7;\n",
" if (u>=-us) curr = 0;\n",
" else curr = Rd*(u+us);\n",
" endif\n",
"endfunction\n",
"function y = circuit(t,u)\n",
" C1 = 1; C2 = 1;\n",
" y = [-1/C1*(Diode(u(1)-10*sin(t))-Diode(u(2)-u(1)));\n",
" 10*cos(t)-1/C2*Diode(u(2)-u(1))];\n",
"endfunction\n",
"t0 = cputime();\n",
"[t,u] = rk45('circuit',0,Tend,u0,1e-5,1e-5); % Runge Kutta adaptiv\n",
"timer = cputime()-t0\n",
"figure(1);\n",
"plot(t,u(:,2)'.')\n",
"grid on; xlabel('time'); ylabel('tension');\n",
"tFix = linspace(0,Tend,100);\n",
"t0 = cputime();\n",
"[tFix,uFix] = ode_Runge('circuit',tFix,u0,1); % Runge Kutta\n",
"timer = cputime()-t0\n",
"figure(2);\n",
"plot(tFix,uFix(:,2),t,u(:,2))\n",
"grid on; xlabel('time'); ylabel('tension');\n",
"legend('u fix','u adapt')"
]
},
{
"cell_type": "code",
"execution_count": 72,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"image/svg+xml": [
""
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"y = -1:0.1:1; %% choose the y values\n",
"v = -1:0.1:1; %% choose the v-values\n",
"\n",
"%% define the function\n",
"function dy = F1(y,v)\n",
" dy = v;\n",
"endfunction\n",
"\n",
"function dv = F2(y,v)\n",
" k = 1; alpha = 0.1;\n",
" dv = -k*y-alpha*v;\n",
"endfunction\n",
"\n",
"%% create the zero vector field\n",
"ny = length(y); nv = length(v);\n",
"V1 = zeros(ny,nv); V2 = zeros(ny,nv);\n",
"\n",
"%% compute the vector field\n",
"for i = 1:ny\n",
" for j = 1:nv\n",
" V1(i,j) = F1(y(i),v(j)); V2(i,j) = F2(y(i),v(j));\n",
" endfor\n",
"endfor\n",
"\n",
"%% choose the scale factor\n",
"scalefactor = 0.4;\n",
"\n",
"figure(1);\n",
"quiver(y,v,V1',V2',scalefactor)\n",
"grid on\n",
"axis('equal')\n",
"axis([min(y),max(y),min(v),max(v)])\n",
"xlabel('position'); ylabel('velocity');\n",
"\n",
"figure(2);\n",
"Vlength = sqrt(V1.^2+V2.^2);\n",
"V1n = V1./Vlength; V2n = V2./Vlength;\n",
"quiver(y,v,V1n',V2n',scalefactor)\n",
"xlabel('position'); ylabel('velocity');\n",
"grid on\n",
"axis('equal')\n",
"axis([min(y),max(y),min(v),max(v)])\n",
"\n",
"\n",
"function dy = ODEPend(y,t)\n",
" k = 1; alpha = 0.1;\n",
" dy = [y(2); -k*y(1)-alpha*y(2)];\n",
"endfunction\n",
"\n",
"T = 0:0.05:25;\n",
"YV = lsode('ODEPend',[0;0.9],T);\n",
"\n",
"figure(3);\n",
"plot(T,YV(:,1),';position;',T,YV(:,2),';velocity;');\n",
"grid on\n",
"xlabel('time'); ylabel('position and velocity');\n",
"\n",
"figure(4)\n",
"hold on\n",
"quiver(y,v,V1n',V2n',scalefactor)\n",
"grid on\n",
"axis('equal')\n",
"axis([min(y),max(y),min(v),max(v)])\n",
"plot(YV(:,1),YV(:,2),'r')\n",
"xlabel('position'); ylabel('velocity');\n",
"hold off"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"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.0"
}
},
"nbformat": 4,
"nbformat_minor": 0
}