水中線状構造物:閉ループ系

\displaystyle{\bf Derivation \#6}
\displaystyle{ \left[\begin{array}{cc|cc} M_{qq}+M_{qq}^a & M_{qp} & 0 & 0\\ M_{pq} & M_{pp}+M_{pp}^b & 0 & 0\\\hline -c_0 I_N  & 0 & I_N  & 0\\ 0 & -d_0 I_N   & 0 & I_N \end{array}\right] \left[\begin{array}{c} \ddot{q} \\ \ddot{p} \\\hline \ddot{q}_a \\ \ddot{p}_b \end{array}\right]}
\displaystyle{+\left[\begin{array}{cc|cc} D_{qq}+D_{qq}^a & D_{qp} & 0 & 0\\ D_{pq} & D_{pp}+D_{qq}^b & 0 & 0\\\hline 0 & 0 & {\cal D}_{qq}^a & 0 \\ 0 & 0 & 0 & {\cal D}_{pp}^b \end{array}\right] \left[\begin{array}{c} \dot{q} \\ \dot{p} \\\hline \dot{q}_a \\ \dot{p}_b \end{array}\right]}
\displaystyle{ +\left[\begin{array}{cc|cc} K_{qq} & 0 & -K_{q_aq_a} & -K_{q_ap_b} \\ 0 & K_{pp} & -K_{p_bq_a} & -K_{p_bp_b} \\\hline 0 & 0 & c_2 I_N  & 0\\ 0 & 0 & 0 & d_2 I_N \end{array}\right] \left[ \begin{array}{c} q \\ p \\\hline q_a \\ p_b \end{array} \right] = \left[\begin{array}{c} f_1 \\ 0 \\\hline 0\\ 0 \end{array}\right] }

\displaystyle{M_{qq}+M_{qq}^a=[m^y_{ij} + \displaystyle\sum_{l,k=1}^NF^y_{ilkj}q_lq_k]_{i,j=1,\cdots,N}+\chi_2I_N}
\displaystyle{M_{qp}=[\displaystyle\sum_{l,k=1}^N N^y_{ilkj} q_l p_k]_{i,j=1,\cdots,N}}
\displaystyle{M_{pq}=[\displaystyle\sum_{l,k=1}^N N^z_{ilkj} p_l q_k]_{i,j=1,\cdots,N}}
\displaystyle{M_{pp}+M_{qq}^b=[m^z_{ij} + \displaystyle\sum_{l,k=1}^NF^z_{ilkj}p_lp_k]_{i,j=1,\cdots,N}+\chi_2I_N}

for i=1:N, for j=1:N
Mqq(i,j)=m_y(i,j);
Mqp(i,j)=0;
Mpq(i,j)=0;
Mpp(i,j)=m_z(i,j);
for l=1:N, for k=1:N
Mqq(i,j)=Mqq(i,j)+F_y(i,l,k,j)*q(l)*q(k);
Mqp(i,j)=Mqp(i,j)+N_y(i,l,k,j)*q(l)*p(k);
Mpq(i,j)=Mpq(i,j)+N_z(i,l,k,j)*p(l)*q(k);
Mpp(i,j)=Mpp(i,j)+F_z(i,l,k,j)*p(l)*p(k);
end, end, end, end
MM=eye(4*N,4*N);
MM(1:2*N,1:2*N)=[Mqq+chi2*eye(N,N) Mqp;Mpq Mpp+chi2*eye(N,N)];
MM(2*N+1:3*N,1:N)=-c_0*eye(N,N);
MM(3*N+1:4*N,N+1:2*N)=-d_0*eye(N,N);

\displaystyle{D_{qq}+D_{qq}^a=[c^y_{ij} +\chi_3\sum_{j=1}^N\int_0^1|v-\dot{\eta}|\phi_i\phi_jd\xi + \displaystyle\sum_{l,k=1}^ND^y_{ilkj}q_lq_k+ \sum_{l,k=1}^NE^y_{ilkj}q_l\dot{q}_k]_{i,j=1,\cdots,N}}
\displaystyle{D_{qp}=[\displaystyle\sum_{l,k=1}^NL^y_{ilkj}q_lp_k+\displaystyle\sum_{l,k=1}^NM^y_{ilkj}q_l\dot{p}_k]_{i,j=1,\cdots,N}}
\displaystyle{D_{pq}=[\displaystyle\sum_{l,k=1}^NL^z_{ilkj}p_lq_k+\displaystyle\sum_{l,k=1}^NM^z_{ilkj}p_l\dot{q}_k]_{i,j=1,\cdots,N} }
\displaystyle{D_{pp}+D_{qq}^b=[c^z_{ij} +\chi_3\sum_{j=1}^N\int_0^1|\dot{\eta}|\psi_i\psi_jd\xi + \displaystyle\sum_{l,k=1}^ND^z_{ilkj}p_lp_k+ \displaystyle\sum_{l,k=1}^NE^z_{ilkj}p_l\dot{p}_k]_{i,j=1,\cdots,N} }
\displaystyle{{\cal D}_{qq}^a=[c_1\displaystyle\sum_{l,k=1}^NS^y_{ijkl}q_{ak}q_{al}]_{i,j=1,\cdots,N}-c_1I_N}
\displaystyle{{\cal D}_{pp}^b=[d_1\displaystyle\sum_{l,k=1}^NS^z_{ijkl}p_{bk}p_{bl}]_{i,j=1,\cdots,N}-d_1I_N}

for i=1:N, for j=1:N
Dqq(i,j)=c_y(i,j)+chi3*c_simp*(abs(V-dy).*f_y(:,i).*f_y(:,j));
Dqp(i,j)=0;
Dpq(i,j)=0;
Dpp(i,j)=c_z(i,j)+chi3*c_simp*(abs(dy).*f_z(:,i).*f_z(:,j));
Dqqa(i,j)=0;
Dppb(i,j)=0;
for l=1:N, for k=1:N
Dqq(i,j)=Dqq(i,j)+D_y(i,l,k,j)*q(l)*q(k)+E_y(i,l,k,j)*q(l)*dq(k);
Dqp(i,j)=Dqp(i,j)+L_y(i,l,k,j)*q(l)*p(k)+M_y(i,l,k,j)*q(l)*dp(k);
Dpq(i,j)=Dpq(i,j)+L_z(i,l,k,j)*p(l)*q(k)+M_z(i,l,k,j)*p(l)*dq(k);
Dpp(i,j)=Dpp(i,j)+D_z(i,l,k,j)*p(l)*p(k)+E_z(i,l,k,j)*p(l)*dp(k);
Dqqa(i,j)=Dqqa(i,j)+c_1*S_y(i,j,k,l)*qa(k)*qa(l);
Dppb(i,j)=Dppb(i,j)+d_1*S_z(i,j,k,l)*pb(k)*pb(l);
end, end, end, end
DD=eye(4*N,4*N);
DD(1:2*N,1:2*N)=[Dqq Dqp;Dpq Dpp];
DD(2*N+1:3*N,2*N+1:3*N)=Dqqa-c_1*eye(N,N);
DD(3*N+1:4*N,3*N+1:4*N)=Dppb-d_1*eye(N,N);

\displaystyle{K_{qq}=[k^y_{ij} + \displaystyle\sum_{l,k=1}^NB^y_{ilkj}q_kq_l+ \displaystyle\sum_{l,k=1}^NH^y_{ilkj}p_kp_l]_{i,j=1,\cdots,N}}
\displaystyle{K_{pp}=[k^z_{ij} + \displaystyle\sum_{l,k=1}^NB^z_{ilkj}p_kp_l+ \displaystyle\sum_{l,k=1}^NH^z_{ilkj}q_kq_l]_{i,j=1,\cdots,N}}
\displaystyle{K_{q_aq_a}=[\chi_4 \displaystyle\sum_{l,k=1}^NS^y_{ijkl}q_{vk}q_{vl}]_{i,j=1,\cdots,N}}
\displaystyle{K_{q_ap_b}=[-\chi_5 \displaystyle\sum_{l,k=1}^NS^y_{ijkl}q_{vk}\dot{q}_{l}]_{i,j=1,\cdots,N}}
\displaystyle{K_{p_bq_a}=[\chi_4 \displaystyle\sum_{l,k=1}^NS^z_{ijkl}p_{vk}p_{vl}]_{i,j=1,\cdots,N}}
\displaystyle{K_{p_bp_b}=[\chi_5 \displaystyle\sum_{l,k=1}^NS^z_{ijkl}p_{vk}\dot{p}_{l}]_{i,j=1,\cdots,N}}

for i=1:N, for j=1:N
Kqq(i,j)=k_y(i,j);
Kpp(i,j)=k_z(i,j);
Kqqa(i,j)=0;
Kqpa(i,j)=0;
Kpqb(i,j)=0;
Kppb(i,j)=0;
for l=1:N, for k=1:N
Kqq(i,j)=Kqq(i,j)+B_y(i,l,k,j)*q(l)*q(k)+H_y(i,l,k,j)*p(l)*p(k);
Kpp(i,j)=Kpp(i,j)+B_z(i,l,k,j)*p(l)*p(k)+H_z(i,l,k,j)*q(l)*q(k);
Kqqa(i,j)=Kqqa(i,j)+chi4*S_y(i,j,k,l)*q(k)*q(l);
Kqpa(i,j)=Kqpa(i,j)-chi5*S_y(i,j,k,l)*q(k)*dp(l);
Kpqb(i,j)=Kpqb(i,j)+chi4*S_z(i,j,k,l)*p(k)*p(l);
Kppb(i,j)=Kppb(i,j)+chi5*S_z(i,j,k,l)*p(k)*dp(l);
end, end, end, end
KK=eye(4*N,4*N);
KK(1:N,1:N)=Kqq;
KK(N+1:2*N,N+1:2*N)=Kpp;
KK(1:2*N,2*N+1:4*N)=-[Kqqa Kqpa;Kpqb Kppb];
KK(2*N+1:3*N,2*N+1:3*N)=c_2*eye(N,N);
KK(3*N+1:4*N,3*N+1:4*N)=d_2*eye(N,N);

\displaystyle{f_1=[\chi_3\int_0^1|v-\dot{\eta}|\phi_iv_id\xi]_{i=1,\cdots,N}}

f1=zeros(N,1);
for i=1:N
f1(i)=chi3*c_simp*(abs(V-dy).*f_y(:,i).*V);
end

●for i,j,k,l=1,\cdots,N

\displaystyle{m^y_{ij}=\int_0^1\phi_i\phi_jd\xi+ \Gamma_m\phi_i(1)\phi_j(1) }
\displaystyle{m^z_{ij}=\int_0^1\psi_i\psi_jd\xi+ \Gamma_m\psi_i(1)\psi_j(1) }

for i=1:N, for j=1:N
m_y(i,j)=c_simp*(f_y(:,i).*f_y(:,j))+G_m*f_y(n_d+1,i)*df_y(n_d+1,j);
end, end
m_z=m_y;

\displaystyle{E^y_{ijk\ell}=\int_0^1 \phi_i\phi'_j\int_0^\xi \phi'_k \phi'_\ell d\xi d\xi +\Gamma_m \phi_i(1)\phi'_j(1) \int_0^1 \phi'_k \phi'_\ell d\xi }
\displaystyle{-\int_0^1 \phi_i\phi''_j\int_\xi^1\int_0^\xi \phi'_k \phi'_\ell d\xi d\xi -\Gamma_m  \int_0^1\phi_i\phi''_j\int_0^1 \phi'_k \phi'_\ell d\xi d\xi}
\displaystyle{E^z_{ijk\ell}=\int_0^1 \psi_i\psi'_j\int_0^\xi \psi'_k \psi'_\ell d\xi d\xi +\Gamma_m \psi_i(1)\psi'_j(1) \int_0^1 \psi'_k \psi'_\ell d\xi }
\displaystyle{-\int_0^1 \psi_i\psi''_j\int_\xi^1\int_0^\xi \psi'_k \psi'_\ell d\xi d\xi -\Gamma_m  \int_0^1\psi_i\psi''_j\int_0^1 \psi'_k \psi'_\ell d\xi d\xi}

for i=1:N, for j=1:N, for k=1:N, for l=1:N
E_y(i,j,k,l)=c_simp*(f_y(:,i).*df_y(:,j).*(simp_u*(df_y(:,k).*df_y(:,l))))…
+G_m*f_y(n_d+1,i)*df_y(n_d+1,j)*c_simp*(df_y(:,k).*df_y(:,l))…
-c_simp*(f_y(:,i).*ddf_y(:,j).*(simp_b*(simp_u*(df_y(:,k).*df_y(:,l)))))…
-G_m*c_simp*(f_y(:,i).*ddf_y(:,j).*(simp_u*(df_y(:,k).*df_y(:,l))));
end, end, end, end
E_z=E_y;

\displaystyle{F^y_{ijk\ell}=\int_0^1 \phi_i\phi'_j\int_0^\xi \phi'_k \phi'_\ell d\xi d\xi +\Gamma_m \phi_i(1)\phi'_j(1) \int_0^1 \phi'_k \phi'_\ell d\xi }
\displaystyle{-\int_0^1 \phi_i\phi''_j\int_\xi^1\int_0^\xi \phi'_k \phi'_\ell d\xi d\xi -\Gamma_m  \int_0^1\phi_i\phi''_j\int_0^1 \phi'_k \phi'_\ell d\xi d\xi}
\displaystyle{F^z_{ijk\ell}=\int_0^1 \psi_i\psi'_j\int_0^\xi \psi'_k \psi'_\ell d\xi d\xi +\Gamma_m \psi_i(1)\psi'_j(1) \int_0^1 \psi'_k \psi'_\ell d\xi }
\displaystyle{-\int_0^1 \psi_i\psi''_j\int_\xi^1\int_0^\xi \psi'_k \psi'_\ell d\xi d\xi -\Gamma_m  \int_0^1\psi_i\psi''_j\int_0^1 \psi'_k \psi'_\ell d\xi d\xi}

for i=1:N, for j=1:N, for k=1:N, for l=1:N
F_y(i,j,k,l)=E_y(i,j,k,l);
end, end, end, end
F_z=F_y;

\displaystyle{M^y_{ijk\ell}= \int_0^1 \phi_i\phi'_j\int_0^\xi \psi'_k \psi'_\ell d\xi d\xi  +\Gamma_m \phi_i(1)\phi'_j(1) \int_0^1 \psi'_k \psi'_\ell d\xi }
\displaystyle{-\int_0^1 \phi_i\phi''_j\int_\xi^1\int_0^\xi \psi'_k \psi'_\ell d\xi d\xi -\Gamma_m  \int_0^1\phi_i\phi''_j\int_0^1 \psi'_k \psi'_\ell d\xi d\xi }
\displaystyle{M^z_{ijk\ell}= \int_0^1 \psi_i\psi'_j\int_0^\xi \phi'_k \phi'_\ell d\xi d\xi  +\Gamma_m \psi_i(1)\psi'_j(1) \int_0^1 \phi'_k \phi'_\ell d\xi }
\displaystyle{-\int_0^1 \psi_i\psi''_j\int_\xi^1\int_0^\xi \phi'_k \phi'_\ell d\xi d\xi -\Gamma_m  \int_0^1\psi_i\psi''_j\int_0^1 \phi'_k \phi'_\ell d\xi d\xi }

for i=1:N, for j=1:N, for k=1:N, for l=1:N
M_y(i,j,k,l)=c_simp*(f_y(:,i).*df_y(:,j).*(simp_u*(df_z(:,k).*df_z(:,l))))…
+G_m*f_y(n_d+1,i)*df_y(n_d+1,j)*c_simp*(df_z(:,k).*df_z(:,l))…
-c_simp*(f_y(:,i).*ddf_y(:,j).*(simp_b*(simp_u*(df_z(:,k).*df_z(:,l)))))…
-G_m*c_simp*(f_y(:,i).*ddf_y(:,j).*(simp_u*(df_z(:,k).*df_z(:,l))));
end, end, end, end
M_z=M_y;

\displaystyle{N^y_{ijk\ell}= \int_0^1 \phi_i\phi'_j\int_0^\xi \psi'_k \psi'_\ell d\xi d\xi  +\Gamma_m \phi_i(1)\phi'_j(1) \int_0^1 \psi'_k \psi'_\ell d\xi }
\displaystyle{-\int_0^1 \phi_i\phi''_j\int_\xi^1\int_0^\xi \psi'_k \psi'_\ell d\xi d\xi -\Gamma_m  \int_0^1\phi_i\phi''_j\int_0^1 \psi'_k \psi'_\ell d\xi d\xi }
\displaystyle{+\int_0^1  (\phi_i \phi_j''\psi_k''\psi_\ell''+\phi_i \phi_j''\psi_k'\psi_\ell'''+3\phi_i \phi_j'\psi_k''\psi_\ell'''+\phi_i \phi_j'^T\psi_k'\psi_\ell'''' )d\xi}
\displaystyle{N^z_{ijk\ell}= \int_0^1 \psi_i\psi'_j\int_0^\xi \phi'_k \phi'_\ell d\xi d\xi  +\Gamma_m \psi_i(1)\psi'_j(1) \int_0^1 \phi'_k \phi'_\ell d\xi }
\displaystyle{-\int_0^1 \psi_i\psi''_j\int_\xi^1\int_0^\xi \phi'_k \phi'_\ell d\xi d\xi -\Gamma_m  \int_0^1\psi_i\psi''_j\int_0^1 \phi'_k \phi'_\ell d\xi d\xi }
\displaystyle{+\int_0^1  (\psi_i \psi_j''\phi_k''\phi_\ell''+\psi_i \psi_j''\phi_k'\phi_\ell'''+3\psi_i \psi_j'\phi_k''\phi_\ell'''+\psi_i \psi_j'\phi_k'\phi_\ell'''' )d\xi}

for i=1:N, for j=1:N, for k=1:N, for l=1:N
N_y(i,j,k,l)=c_simp*(f_y(:,i).*df_y(:,j).*(simp_u*(df_z(:,k).*df_z(:,l))))…
+G_m*f_y(n_d+1,i)*df_y(n_d+1,j)*c_simp*(df_z(:,k).*df_z(:,l))…
-c_simp*(f_y(:,i).*ddf_y(:,j).*(simp_b*(simp_u*(df_z(:,k).*df_z(:,l)))))…
-G_m*c_simp*(f_y(:,i).*ddf_y(:,j).*(simp_u*(df_z(:,k).*df_z(:,l))))…
+c_simp*(f_y(:,i).*ddf_y(:,j).*ddf_z(:,k).*ddf_z(:,l)…
+f_y(:,i).*ddf_y(:,j).*df_z(:,k).*dddf_z(:,l)…
+3*f_y(:,i).*df_y(:,j) .*ddf_z(:,k).*dddf_z(:,l)…
+f_y(:,i).*df_y(:,j) .*df_y(:,k).*ddddf_y(:,l));
end, end, end, end
N_z=N_y;

\displaystyle{S^y_{ijkl}=\int_0^1 \phi_i\phi_j\phi_k\phi_l d\xi}
\displaystyle{S^z_{ijkl}=\int_0^1 \psi_i\psi_j\psi_k\psi_l d\xi}

for i=1:N, for j=1:N, for k=1:N, for l=1:N
S_y(i,j,k,l)=c_simp*(f_y(:,i).*f_y(:,j).*f_y(:,k).*f_y(:,l));
end, end, end, end
S_z=S_y;

\displaystyle{c^y_{ij}= \int_0^1 2 \sqrt{\beta}{u} \phi_i \phi'_j d\xi }
\displaystyle{c^z_{ij}= \int_0^1 2 \sqrt{\beta}{u} \psi_i \psi'_j d\xi }

for i=1:N, for j=1:N
c_y(i,j)=2*sqrt(beta)*c_simp*(U.*f_y(:,i).*df_y(:,j));
end, end
c_z=c_y;

\displaystyle{k^y_{ij}= \int_0^1\gamma\phi_i\phi'_jd\xi+\gamma\Gamma_p\phi_i(1)\phi'_j(1) }
\displaystyle{- \int_0^1 (\gamma(1-\xi+\Gamma_p)} \displaystyle{-\sqrt{\beta}{\dot{u}}(1-\xi) } \displaystyle{) \phi_i \phi_j''d\xi } \displaystyle{+\int_0^1 {u^2} \phi_i\phi''_j  d\xi } \displaystyle{+ \int_0^1 \phi_i \phi''''_j  d\xi}
\displaystyle{k^z_{ij}= \int_0^1\gamma\psi_i\psi'_jd\xi+\gamma\Gamma_p\psi_i(1)\psi'_j(1) }
\displaystyle{- \int_0^1 (\gamma(1-\xi+\Gamma_p)} \displaystyle{-\sqrt{\beta}{\dot{u}}(1-\xi) } \displaystyle{) \psi_i \psi_j''d\xi } \displaystyle{+\int_0^1  {u^2} \psi_i\psi''_j  d\xi } \displaystyle{+ \int_0^1 \psi_i \psi''''_j  d\xi}

for i=1:N, for j=1:N
k_y_0(i,j)=c_simp*(gamma*f_y(:,i).*df_y(:,j))+gamma*G_p*f_y(n_d+1,i).*df_y(n_d+1,j)…
-c_simp*(gamma*(1-xi+G_p))+c_simp*(f_y(:,i).*ddddf_y(:,j));
end, end

for i=1:N, for j=1:N
k_y(i,j)= k_y_0(i,j)+sqrt(beta)*c_simp*(dU.*(1-xi).*f_y(:,i).*ddf_y(:,j))…
+c_simp*(U.^2.*f_y(:,i).*ddf_y(:,j));
end, end
k_z=k_y;

\displaystyle{B^y_{ijk\ell}={1\over 2}\gamma \int_0^1 \phi_i\phi'_j\phi'_k \phi'_\ell d\xi +{1\over 2}\gamma\Gamma_p \phi_i(1)\phi'_j(1) \phi'_k(1) \phi'_\ell(1)}
\displaystyle{-{3\over 2} \int_0^1 (\gamma(1-\xi+\Gamma_p)} \displaystyle{-\sqrt{\beta}{\dot{u}}(1-\xi) } \displaystyle{) \phi_i \phi_j'' \phi'_k \phi'_\ell d\xi }
\displaystyle{+\int_0^1 \sqrt{\beta}{\dot{u}}\phi_i \phi_j'' \int_\xi^1  \phi'_k \phi'_\ell d\xi d\xi +\int_0^1  {u^2} ( \phi_i\phi_j' \phi'_k \phi''_\ell-  \phi_i\phi''_j\int_\xi^1 \phi'_k \phi''_\ell d\xi) d\xi }
\displaystyle{+\int_0^1  (\phi_i \phi_j''\phi_k''\phi_\ell''+4\phi_i \phi_j'\phi_k''\phi_\ell'''+\phi_i \phi_j''\phi_k''\phi_\ell'''' )d\xi}
\displaystyle{B^z_{ijk\ell}={1\over 2}\gamma \int_0^1 \psi_i\psi'_j\psi'_k \psi'_\ell d\xi +{1\over 2}\gamma\Gamma_p \psi_i(1)\psi'_j(1) \psi'_k(1) \psi'_\ell(1)}
\displaystyle{-{3\over 2} \int_0^1 (\gamma(1-\xi+\Gamma_p)} \displaystyle{-\sqrt{\beta}{\dot{u}}(1-\xi) } \displaystyle{) \psi_i \psi_j'' \psi'_k \psi'_\ell d\xi }
\displaystyle{+\int_0^1 \sqrt{\beta}{\dot{u}}\psi_i \psi_j'' \int_\xi^1  \psi'_k \psi'_\ell d\xi d\xi +\int_0^1  {u^2} ( \psi_i\psi_j' \psi'_k \psi''_\ell-  \psi_i\psi''_j\int_\xi^1 \psi'_k \psi''_\ell d\xi) d\xi }
\displaystyle{+\int_0^1  (\psi_i \psi_j''\psi_k''\psi_\ell''+4\psi_i \psi_j'\psi_k''\psi_\ell'''+\psi_i \psi_j''\psi_k''\psi_\ell'''' )d\xi}

for i=1:N, for j=1:N, for k=1:N, for l=1:N
B_y_0(i,j,k,l)=1/2*gamma*c_simp*(f_y(:,i).*df_y(:,j).*df_y(:,k).*df_y(:,l))…
+1/2*gamma*G_p*f_y(n_d+1,i).*df_y(n_d+1,j).*df_y(n_d+1,k).*df_y(n_d+1,l)…
-3/2*c_simp*(gamma*(1-xi+G_p))…
+c_simp*(f_y(:,i).*ddf_y(:,j).*ddf_y(:,k).*ddf_y(:,l)…
+4*f_y(:,i).*df_y(:,j) .*ddf_y(:,k).*dddf_y(:,l)…
+f_y(:,i).*ddf_y(:,j).*ddf_y(:,k).*ddddf_y(:,l));
end, end, end, end

for i=1:N, for j=1:N, for k=1:N, for l=1:N
B_y(i,j,k,l)= B_y_0(i,j,k,l)…
+3/2*sqrt(beta)*c_simp*(dU.*(1-xi).*f_y(:,i).*ddf_y(:,j).*df_y(:,k).*df_y(:,l))…
+sqrt(beta)*c_simp*(dU.*f_y(:,i).*ddf_y(:,j).*(simp_b*(df_y(:,k).*df_y(:,l))))…
+c_simp*(U.^2.*(f_y(:,i).*df_y(:,j).*df_y(:,k).*ddf_y(:,l)…
-f_y(:,i).*ddf_y(:,j).*(simp_b*(df_y(:,k).*ddf_y(:,l)))));
end, end, end, end
B_z=B_y;

\displaystyle{H^y_{ijk\ell}={1\over 2}\gamma \int_0^1 \phi_i\phi'_j\psi'_k \psi'_\ell d\xi +{1\over 2}\gamma\Gamma_p \phi_i(1)\phi'_j(1) \psi'_k(1) \psi'_\ell(1)}
\displaystyle{-{1\over 2} \int_0^1(\gamma(1-\xi+\Gamma_p)} \displaystyle{-\sqrt{\beta}{\dot{u}}(1-\xi) } \displaystyle{) \phi_i \phi_j'' \psi'_k \psi'_\ell d\xi  }
\displaystyle{- \int_0^1 (\gamma(1-\xi+\Gamma_p)} \displaystyle{-\sqrt{\beta}{\dot{u}}(1-\xi) } \displaystyle{) \phi_i \phi_j' \psi'_k \psi''_\ell d\xi  }
\displaystyle{-\int_0^1 \sqrt{\beta}{\dot{u}}\phi_i \phi_j'' \int_\xi^1  \psi'_k \psi'_\ell d\xi d\xi +\int_0^1  {u^2} ( \phi_i\phi_j' \psi'_k \psi''_\ell -  \phi_i\phi''_j\int_\xi^1 \psi'_k \psi''_\ell d\xi) d\xi }
\displaystyle{H^z_{ijk\ell}={1\over 2}\gamma \int_0^1 \psi_i\psi'_j\phi'_k \phi'_\ell d\xi +{1\over 2}\gamma\Gamma_p \psi_i(1)\psi'_j(1) \phi'_k(1) \phi'_\ell(1)}
\displaystyle{-{1\over 2} \int_0^1(\gamma(1-\xi+\Gamma_p)} \displaystyle{-\sqrt{\beta}{\dot{u}}(1-\xi) } \displaystyle{) \psi_i \psi_j'' \phi'_k \phi'_\ell d\xi  }
\displaystyle{- \int_0^1 (\gamma(1-\xi+\Gamma_p)} \displaystyle{-\sqrt{\beta}{\dot{u}}(1-\xi) } \displaystyle{) \psi_i \psi_j' \phi'_k \phi''_\ell d\xi  }
\displaystyle{-\int_0^1 \sqrt{\beta}{\dot{u}}\psi_i \psi_j'' \int_\xi^1  \phi'_k \phi'_\ell d\xi d\xi +\int_0^1  {u^2} ( \psi_i\psi_j' \phi'_k \phi''_\ell -  \psi_i\psi''_j\int_\xi^1 \phi'_k \phi''_\ell d\xi) d\xi }

for i=1:N, for j=1:N, for k=1:N, for l=1:N
H_y_0(i,j,k,l)=1/2*gamma*c_simp*(f_y(:,i).*df_y(:,j).*df_y(:,k).*df_y(:,l))…
+1/2*gamma*G_p*f_y(n_d+1,i).*df_y(n_d+1,j).*df_y(n_d+1,k).*df_y(n_d+1,l)…
-1/2*c_simp*(gamma*(1-xi+G_p).*f_y(:,i).*ddf_y(:,j).*df_y(:,k).*df_y(:,l))…
-c_simp*(gamma*(1-xi+G_p).*f_y(:,i).*df_y(:,j).*df_z(:,k).*ddf_z(:,l));
end, end, end, end

for i=1:N, for j=1:N, for k=1:N, for l=1:N
H_y(i,j,k,l)= H_y_0(i,j,k,l)…
+1/2*sqrt(beta)*c_simp*(dU.*(1-xi).*f_y(:,i).*ddf_y(:,j).*df_y(:,k).*df_y(:,l)…
+sqrt(beta)*(1-xi).*f_y(:,i).*df_y(:,j).*df_z(:,k).*ddf_z(:,l)…
+sqrt(beta)*f_y(:,i).*ddf_y(:,j).*(simp_b*(df_z(:,k).*df_z(:,l))))…
+c_simp*(U.^2.*(f_y(:,i).*df_y(:,j).*df_z(:,k).*ddf_z(:,l)…
-f_y(:,i).*ddf_y(:,j).*(simp_b*(df_y(:,k).*ddf_y(:,l)))));
end, end, end, end
H_z=H_y;

\displaystyle{D^y_{ijk\ell}=\int_0^1  2 \sqrt{\beta}{u} ( \phi_i\phi_j' \phi'_k \phi'_\ell-  \phi_i\phi''_j\int_\xi^1  \phi'_k \phi'_\ell d\xi) d\xi}
\displaystyle{D^z_{ijk\ell}=\int_0^1  2 \sqrt{\beta}{u} ( \psi_i\psi_j' \psi'_k \psi'_\ell-  \psi_i\psi''_j\int_\xi^1  \psi'_k \psi'_\ell d\xi) d\xi}

for i=1:N, for j=1:N, for k=1:N, for l=1:N
D_y(i,j,k,l)=2*sqrt(beta)*c_simp*(U.*(f_y(:,i).*df_y(:,j).*df_y(:,k).*df_y(:,l)…
-f_y(:,i).*ddf_y(:,j).*(simp_b*(df_y(:,k).*df_y(:,l)))));
end, end, end, end
D_z=D_y;

\displaystyle{L^y_{ijk\ell}=\int_0^1  2 \sqrt{\beta}{u} ( \phi_i\phi_j'  \psi'_k \psi'_\ell-  \phi_i\phi''_j\int_\xi^1  \psi'_k \psi'_\ell d\xi) d\xi}
\displaystyle{L^z_{ijk\ell}=\int_0^1  2 \sqrt{\beta}{u} ( \psi_i\psi_j'  \phi'_k \phi'_\ell-  \psi_i\psi''_j\int_\xi^1  \phi'_k \phi'_\ell d\xi) d\xi}

for i=1:N, for j=1:N, for k=1:N, for l=1:N
L_y(i,j,k,l)=2*sqrt(beta)*c_simp*(U.*(f_y(:,i).*df_y(:,j).*df_z(:,k).*df_z(:,l)…
-f_y(:,i).*ddf_y(:,j).*(simp_b*(df_z(:,k).*df_z(:,l)))));
end, end, end, end
L_z=L_y;