Cobweb plot

While conducting multi variable analysis, it may be useful to see how the output changes due to change in inputs. For example humidity, temperature and water-level may change the deflection of dam. Cobweb plot combined with normal distribution will help to visualize the risks easily. An example of such plot is shown below.

An octave (or matlab) code to generate such plot is given below. Basically, it requires to plot each entity separatly. The data should be normalized to set equal scale.

%shaded plot
clear all;clf; % clear figure (just in case something was there before)

epdata=[-0.0008959	-0.00059726	-0.0013775	-0.00080342	-0.002345	-0.0012454	-0.0014801	-0.00097711	-0.00096319	-0.00030882	-0.0015608	-0.00099932	-0.00082923	-0.00105	-0.0013915	-0.00021562	-0.0011739	-0.0015718	-0.00087249	-0.0016794	-0.00036352	-0.00059309	-0.0011483	-0.0010326	-0.00034531	-0.0020304	-0.0014554	-0.00092798	-0.0012557	-0.0016694	-0.0010306	-0.001245	-0.0013191	-0.0016508	-0.0010848	-0.0018287	-0.0010912	-0.0014441	-0.0012839	-0.00089405	-0.0024034	-0.0015758	-0.00026981	-0.0014717	-0.0023293	-0.0020914	-0.00076314	-0.00056843	-0.0014572	-0.00083516	-0.0012168	-0.0013884	-0.0012994	-0.0010078	-0.00045422	-0.0015285	-0.00061293	-0.00091066	-0.0014499	-0.0010975	-0.001501	-0.0012706	-0.0019573	-0.0019029	-0.0013711	-0.0019286	-0.00099036	-0.0016533	-0.001049	-0.0015944	-0.0016858	-0.0014946	0.000314	-0.0016982	-0.0012176	-0.0016377	-0.0012109	-0.0015727	-0.001435	-0.0010661	-0.0012172	-0.0020125]
tcdata=[31.007	17.279	71.559	42.942	30.757	64.937	56.682	34.429	48.118	51.053	46.382	70.517	49.346	56.116	38.708	39.931	61.857	52.765	58.071	55.669	34.416	38.345	36.623	42.745	66.059	83.693	40.504	62.218	64.917	20.47	65.233	37.162	69.894	27.246	3.5394	65.416	65.68	43.018	95.189	59.176	59.944	56.573	53.296	38.949	41.054	39.729	28.056	43.599	59.708	46.291	49.449	48.566	57.104	51.747	68.345	48.032	63.941	41.263	26.673	48.541	31.849	56.63	39.846	49.557	80.801	37.334	37.459	36.248	57.355	39.958	32.182	24.017	75.422	61.597	43.402	11.902	58.988	42.176	17.418	44.48	56.916	41.854]
tldata=[61.235	77.001	86.213	75.22	43.764	67.211	76.459	52.866	77.028	54.373	34.241	90.087	80.787	60.298	75.159	107.12	84.327	39.385	104.76	66.38	62.929	42.629	88.112	83.559	59.862	85.192	50.859	55.306	66.787	56.848	82.743	95.426	112.39	44.012	70.265	54.671	76.631	76.642	88.841	39.25	56.877	60.893	104.51	59.836	63.706	71.222	63.624	32.153	70.074	79.584	72.205	72.182	73.233	67.998	80.533	40.931	63.258	42.223	79.662	81.851	52.266	80.383	83.354	57.363	93.174	64.228	87.165	93.264	79.508	67.241	98.564	86.5	37.035	59.61	64.813	68.867	72.937	81.84	61.348	92.968	57.552	64.609]
dpdata=[-0.016616	0.0012712	-0.018897	-0.01016	-0.092141	-0.024361	-0.023951	-0.026684	-0.013089	-0.0056231	-0.065296	-0.011672	-0.0096546	-0.023527	-0.018468	0.0022514	-0.015527	-0.056125	-0.0052227	-0.0348	-0.0046483	-0.020469	-0.0072833	-0.010084	-0.0052599	-0.030438	-0.044018	-0.022032	-0.024795	-0.031457	-0.014002	-0.0055563	-0.0096733	-0.063971	0.003394	-0.046408	-0.017231	-0.019991	-0.016868	-0.028567	-0.061399	-0.036514	0.00094998	-0.034091	-0.050397	-0.03461	-0.011083	-0.023391	-0.027504	-0.0096521	-0.020398	-0.023604	-0.022343	-0.018593	-0.0046534	-0.055475	-0.011241	-0.03279	-0.0084095	-0.013231	-0.043265	-0.018161	-0.020591	-0.04901	-0.016887	-0.038989	-0.0064535	-0.0089957	-0.014825	-0.029465	-0.0045653	-0.0031735	0.013363	-0.040035	-0.024299	-0.00078227	-0.020854	-0.018003	-0.015243	-0.0074379	-0.02939	-0.042169]
% normalize
ep=epdata/min(epdata);
tc=tcdata/max(tcdata);
tl=tldata/max(tldata);
dp=dpdata/min(dpdata);
x=[0 1 2 3.5]*7;
x1=ones(1,length(ep))*x(1);
x2=ones(1,length(ep))*x(2);
x3=ones(1,length(ep))*x(3);
x4=ones(1,length(ep))*x(4);

%calculate mean
m_ep=mean(ep);
m_tc=mean(tc);
m_tl=mean(tl);
m_dp=mean(dp);
%standard deviation
sd_ep=std(ep);
sd_tc=std(tc);
sd_tl=std(tl);
sd_dp=std(dp);

%coords to plot normal distribution
nd_x_ep=[m_ep-4*sd_ep:sd_ep/20:m_ep+4*sd_ep]
nd_y_ep=1/sd_ep/(2*pi())^0.5*exp((nd_x_ep-m_ep).^2/-2/sd_ep^2)

nd_x_tc=[m_tc-4*sd_tc:sd_tc/20:m_tc+4*sd_tc]
nd_y_tc=1/sd_tc/(2*pi())^0.5*exp((nd_x_tc-m_tc).^2/-2/sd_tc^2)

nd_x_tl=[m_tl-4*sd_tl:sd_tl/20:m_tl+4*sd_tl]
nd_y_tl=1/sd_tl/(2*pi())^0.5*exp((nd_x_tl-m_tl).^2/-2/sd_tl^2)

nd_x_dp=[m_dp-4*sd_dp:sd_dp/20:m_dp+4*sd_dp]
nd_y_dp=1/sd_dp/(2*pi())^0.5*exp((nd_x_dp-m_dp).^2/-2/sd_dp^2)

figure 1
% plot vertical lines
hold on;
plot (x1,ep,"o"); 
plot (x2,tc,"o"); 
plot (x3,tl,"o"); 
plot (x4,dp,"o"); 
%plot connecting lines
p=[ep; tc; tl; dp ];
plot (x, p, "color",[.8 .8 .8] );

%plot the normal disribution
patch(nd_y_ep*2+x(1),nd_x_ep,'facecolor',[1 1 0.8],'edgecolor',[1 0 0], 'facealpha',0.3)
patch(nd_y_tc*2+x(2),nd_x_tc,'facecolor',[1 1 0.8],'edgecolor',[1 0 0], 'facealpha',0.3)
patch(nd_y_tl*2+x(3),nd_x_tl,'facecolor',[1 1 0.8],'edgecolor',[1 0 0], 'facealpha',0.3)
patch(-nd_y_dp*2+x(4),nd_x_dp,'facecolor',[1 1 0.8],'edgecolor',[0 0 1], 'facealpha',0.3,"linewidth", 3)

%yaxis lines
ymin=-0.5;ymax=1.5;
plot ([x(1) x(1)],[ymin ymax],"color",[0 0 0] );
plot ([x(2) x(2)],[ymin ymax],"color",[0 0 0] );
plot ([x(3) x(3)],[ymin ymax],"color",[0 0 0] );
plot ([x(4) x(4)],[ymin ymax],"color",[0 0 0] );
%other properties
ylim([ymin ymax])
axis off;
set(gca, "linewidth", 2, "fontsize", 12);
%labels
%xlabels
text(x(1)+0.3,ymin,"Strain","color","black","fontsize", 18)
text(x(2)+0.3,ymin,"TC","color","black","fontsize", 18)
text(x(3)+0.3,ymin,"TL","color","black","fontsize", 18)
text(x(4)-3.5,ymin,"Max. Deflection","color","black","fontsize", 18)
%informations
text(x(1)+0.3,ymax-.2,["Mean=" num2str(-0.0012, "%5.2e") char(10) "SD=" num2str(0.005, "%5.2e")],"color","black","fontsize", 15)
text(x(2)+0.3,ymax-.2,["Mean=" num2str(50, "%5.0f") char(10) "SD=" num2str(15, "%5.0f") " years"],"color","black","fontsize", 15)
text(x(3)+0.3,ymax-.2,["Mean=" num2str(70, "%5.0f") char(10) "SD=" num2str(20, "%5.0f") " years"],"color","black","fontsize", 15)
text(x(4)-3.5,ymax-.2,["Mean=" num2str(mean(dpdata), "%5.2e") char(10) "SD=" num2str(std(dpdata), "%5.2e") " m"],"color","black","fontsize", 15)

hold off;
print -dpng figure1.png

Advertisement

Derivation of Euler’s forumula for imaginary numbers [e(x)=coz(x)+i*sin(x)]

Derivation of Euler’s forumula for imaginary numbers
a+ib=r(cosθ+ i sinθ)

From Tylor’s series:
F(x)=F(0)+x*f’(x)/1!+x^2*f’’(0)/2!+……+x^n*fn(0)/n!+….
Therefore,
Sin(x)=x-x^3/3!+x^5/5!-x^7/7!+….
Cos(x)=1-x^2/2!+x^4/4!-x^6/6!+….
e(x)=1+x+x^2/2!+x^3/3!+..

Now,
E(ix)=1+ix+(ix)^2/2!+(ix)^3/3!+(ix)^4/4!+(ix)^5/5!+…
E(ix)=(1-x^2/2!+x^4/4!-x^6/6!+…)+i(x-x^3/3!+x^5/5!-x^7/7!+…)
=cos(x)-i*sin(x)

the video for derivation is here

for derivation of Taylor’s series:

Derivation for sum of angles

Cos(A+B)=OP/OS=(OQ-PQ)/OS
	=OQ/OS-PQ/OS
	=OQ/OR*OR/OS-PQ/SR*SR/OS
	=CosA*CosB-SinA*SinB …………(since PQ=TR)
Sin(A+B)=SP/OS=(ST+TP)/OS
	=ST/OS+TP/OS
	=ST/SR*SR/OS+TP/OR*OR/OS
	=CosA*SinB +SinA*CosB	………..(since TP=RQ)
Tan(A+B)=Sin(A+B)/Cos(A+B)
	=(SinA*CosB+CosA*SinB)/(CosA*CosB-SinA*SinB)
	=(TanA+TanB)/(1-TanA*TanB)………(dividing both side by CosA*CosB) 

For A-B note that sin is odd function and cos is even function, therefore, Sin(-A)=-SinA and Cos(-A)=Cos(A). Graphically this means the cos function is symmetric over x axis while sine function is not so.

Radius of earth (classical trignometry)

Question: Suppose that, while lying on a beach near the equator watching the sun set over a calm ocean, you start a stopwatch just as the top of the sun disappears. You then stand, elevating your eyes by a height H=1.7m, and stop the watch when the top of the sun again disappears. If the elapsed time is t=11.1 secs, what is the radius r of the earth?

Solution
One day=24hrs=86400sec

The earht rotate 2*pi() in one day=2*pi=6.2832 radians

angle traversed in 11.1sec=0.000807215 radians

observation height=1.7m

R=(h*cosX)/(1-cosX)=5.22E+06

Download calculation file

Solving Laplace equation in MS Excel for two dimensional solution

Laplace equation is second order derivative of the form shown below.

where phi is a potential function.

This equation is used to describe the behavior of electric, gravitational, and fluid potentials. In the study of heat conduction, the Laplace equation is the steady-state heat equation. This equation also describes seepage underneath the dam.

Direct solution of this equation is very tedious( and uninteresting). Numerical technique gives another approach to solve this equation. In numerical method, the solution space is converted in discrete number of points (a two dimensional array). Any point ‘p’ inside the solution matrix is the average of its surrounding points (4 points if 4 direction is considerer, 6 points if diagonal elements are included). The known points will be the potential(Phi value) at the boundary. Iteration has to done until desired accuracy is obtained. (similar to playing Minesweeper)

Excel sheet is nothing but a huge matrix.Thus providing ideal place to play with matrix. To solve Laplace equation, first of all each cell should be given a fixed width, i.e. say each cell will be 0.5m. After that, draw boundary, by putting the known phi value at the boundary. Any cells within the boundary will now given a formula which says =average(all surrounding cells). By doing so you will create cyclic function which will run automatically until desired level of accuracy is met. You must enable iterative calculation and set the desired accuracy in excel in its option menu.

I will demonstrate this with an example of seepage problem underneath the dam(green) as shown in figure. The dam has a pile in the upstream toe. There is seepage treatment 25m u/s and 30m d/s of dam. The base width of dam is 10m. The pile is 7m deep sheet pile (where seepage=0). The u/s water head is 5m and d/s water head is 1m. We are interested in distribution of water head inside the soil ( shown in yellow). As 2 dimensional Laplace equation will describe the flow field in this problem, we will solve it in excel sheet. Let us suppose at 25m below the dam foundation there is impermeable layer (another impermeable boundary).

Each cell is considered to be1.0m wide. Red cells denotes the boundary. The pile is also hydraulic boundary. The green is dam (done to show similarity with picture). Now a formula is written in any cell stating it as average of 4 surrounding cells (for e.g I10=AVERAGE(I9,H10,I11,J10)) and copy the cell in yellow zone. The boundary condition at u/s is phi=5m water head and at d/s is phi=1m of water head. After you copy formula to each of yellow cell, excel starts to calculate. For this iterative calculation has to be turned on (check Excel>options>formula). The higher the number of iteration, more accurate the result. Sometimes you may have to hit “Calculate sheet” button multiple times to get the result. A chart will show equipotential lines. A sample is shown in figure below. The excel sheet can be downloaded here.