Question: The procedure can be operation when MM:=evalm(Hi.EE.EEE.H0.H1.H2.H3.H4):,but it cann't be opration when MM:=evalm(Hi.EE.EEE.H0.H1.H2.H3.H4.H5):,why?

restart:
Digits:=30;

h0:=0.156;
d:=0.32*h0;
l:=1;
h1:=h0-d;
h2:=h0+d;
h3:=0.5*h0;
g:=9.8;
d1:=1;
Term:=5;
Num:=150:
n:=1:
l1:=l/n;
p:=2;

for N from 1 to Num do
lambda:=2*n*Pi/l:## N1 wei sha ba tiao shu 
k0:=evalf(0.5*Pi+2*(N-1)*Pi/(Num-1)):
tau0:=evalf(k0*h0):
omega:=evalf((g*k0*tanh(k0*h0))^(1/2)):
E:=g/(omega)^2:
k1:=abs(fsolve(omega^2=g*k*tanh(k*h1),k)):
tau1:=evalf(k1*h1):
k2:=abs(fsolve(omega^2=g*k*tanh(k*h2),k)):
tau2:=evalf(k2*h2):
k3:=abs(fsolve(omega^2=g*k*tanh(k*h3),k)):
tau3:=evalf(k3*h3):

F:=tau->tanh(tau1)+tau*add((eval(diff(tanh(tau),tau$s),tau=tau1))/s!*(tau-tau1)^(s-1),s=1..10);##F1(K)
T:=tau->tanh(tau2)+tau*add((eval(diff(tanh(tau),tau$s),tau=tau2))/s!*(tau-tau2)^(s-1),s=1..10);##F2(K)

P:=tau->(sinh(2*tau)-2*tau*cosh(2*tau))/(2*tau+sinh(2*tau))^2;##P
Q:=tau->(16*tau^4+32*tau^3*sinh(2*tau)-9*sinh(2*tau)*sinh(4*tau)+12*tau*(tau+sinh(2*tau))*((cosh(2*tau))^2-2*cosh(2*tau)+3))/3/(2*tau+sinh(2*tau))^4; ##Q
A:=unapply(taylor(2*(tau-tau1)/sinh(2*tau)-tanh(tau)*(h0-E*tau*tanh(tau))*(2*tau+sinh(2*tau))/((E*tau*tanh(tau)-h2)*F(tau)*sinh(2*tau)),tau=tau1,Term+1),tau);##A(K)
B:=unapply(taylor((1+2*tau/sinh(2*tau))^2*(-(tau-tau1)/lambda^2/E/(E*tau*tanh(tau)-h2)/F(tau)-(h0-E*tau*tanh(tau))*(tau-tau1)*tanh(tau)/(E*tau*tanh(tau)-h2)/F(tau)*P(tau)+(tau-tau1)^2*Q(tau)),tau=tau1,Term+1),tau);##B(K)


for j from 0 to Term do
a[j]:=coeff(A(tau),tau-tau1,j):
b[j]:=coeff(B(tau),tau-tau1,j):
end do;

for m from 1 to Term do
f1[0]:=1;
f2[0]:=1;
f1[m]:=-(add(f1[m-i]*((m-i)*a[i]+b[i]),i=1..m))/(m*(m-1+1/2));##pm(z1)
f2[m]:=-(add(f2[m-i]*((m-i+1/2)*a[i]+b[i]),i=1..m))/((m+1/2)*(m+1/2-1+1/2));##qm(Z2)
end do;


CC:=unapply(taylor(2*(tau-tau2)/sinh(2*tau)-tanh(tau)*(h0-E*tau*tanh(tau))*(2*tau+sinh(2*tau))/(E*tau*tanh(tau)-h1)/T(tau)/sinh(2*tau),tau=tau2,Term+1),tau);
DD:=unapply(taylor((1+2*tau/sinh(2*tau))^2*(-(tau-tau2)/lambda^2/E/(E*tau*tanh(tau)-h1)/T(tau)-(h0-E*tau*tanh(tau))*(tau-tau2)*tanh(tau)/(E*tau*tanh(tau)-h1)/T(tau)*P(tau)+(tau-tau2)^2*Q(tau)),tau=tau2,Term+1),tau);

j:='j':
for j from 0 to Term do
cc[j]:=coeff(CC(tau),tau-tau2,j):
dd[j]:=coeff(DD(tau),tau-tau2,j):
end do;

i:='i':
f3[0]:=1;
f4[0]:=1;
m:='m':
for m from 1 to Term do
f3[m]:=-(add(f3[m-i]*((m-i)*cc[i]+dd[i]),i=1..m))/(m*(m-1+1/2));
f4[m]:=-(add(f4[m-i]*((m-i+1/2)*cc[i]+dd[i]),i=1..m))/((m+1/2)*(m+1/2-1+1/2));
end do:

xi11:=unapply(add(f1[j1]*(tau-tau1)^(j1),j1=0..m-1),tau);
xi12:=unapply(add(f2[j2]*(tau-tau1)^(j2+1/2),j2=0..m-1),tau);
xi21:=unapply(add(f3[j3]*(tau-tau2)^(j3),j3=0..m-1),tau);
xi22:=unapply(add(f4[j4]*(tau-tau2)^(j4+1/2),j4=0..m-1),tau);


u0:=evalf(g*tanh(tau0)*(1+2*tau0/(sinh(2*tau0)))/(2*k0));
u01:=evalf(g*tanh(tau3)*(1+2*tau3/(sinh(2*tau3)))/(2*k3));
u1:=evalf(g*(1-(tanh(tau0))^2)*(sinh(2*tau0)-2*tau0*cosh(2*tau0))/(4*(2*tau0+sinh(2*tau0))));
H:=evalf((1+2*tau0/sinh(2*tau0))/(-lambda*k0*d));
delta00:=evalf((lambda*d*u1/u0+I*k0)*H);
delta01:=evalf((lambda*d*u1/u0-I*k0)*H);
delta11:=evalf((lambda*d*u1/u0+I*k0)*H*exp(I*k0*l));
delta12:=evalf((lambda*d*u1/u0-I*k0)*H*exp(-I*k0*l));
delta21:=evalf(exp(I*k0*(l1)));
delta22:=evalf(u0*k0*exp(I*k0*(l1)));
delta31:=evalf(exp(-I*k0*(l1)));
delta32:=evalf(-u0*k0*exp(-I*k0*(l1)));
delta41:=evalf(exp(I*k3*(l1)));
delta42:=evalf(u01*k3*exp(I*k3*(l1)));
delta51:=evalf(exp(-I*k3*(l1)));
delta52:=evalf(-u01*k3*exp(-I*k3*(l1)));
delta61:=evalf(exp(I*k3*(l1+d1)));
delta62:=evalf(u01*k3*exp(I*k3*(l1+d1)));
delta71:=evalf(exp(-I*k3*(l1+d1)));
delta72:=evalf(-u01*k3*exp(-I*k3*(l1+d1)));
Y11 := evalf(exp(k0*(l1 + d1)*I));
Y12 := evalf(-exp(-k0*(l1 + d1)*I));
Y21 := evalf(u0*k0*exp(k0*(l1 + d1)*I));
Y22 := evalf(-u0*k0*exp(-k0*(l1 + d1)*I));

G11 := evalf(exp(k0*(2*l1 + d1)*I));
G12 := evalf(-exp(-k0*(2*l1 + d1)*I));
G21 := evalf(u0*k0*exp(k0*(2*l1 + d1)*I));
G22 := evalf(-u0*k0*exp(-k0*(2*l1 + d1)*I));

W11 := evalf(exp(k3*(2*l1 + d1)*I));
W12 := evalf(-exp(-k3*(2*l1 + d1)*I));
W21 := evalf(u01*k3*exp(k3*(2*l1 + d1)*I));
W22 := evalf(-u01*k3*exp(-k3*(2*l1 + d1)*I));

Z11 := evalf(exp(k3*(2*l1 + 2*d1)*I));
Z12 := evalf(-exp(-k3*(2*l1 + 2*d1)*I));
Z21 := evalf(u01*k3*exp(k3*(2*l1 + 2*d1)*I));
Z22 := evalf(-u01*k3*exp(-k3*(2*l1 + 2*d1)*I));


V11:=evalf(exp(I*k0*2*(l1+d1)));
V21:=evalf(u0*k0*exp(I*k0*2*(l1+d1)));

delta91:=evalf(exp(I*k0*l));
delta92:=evalf(exp(-I*k0*l));


Hi:=(Matrix([[1,1],[delta00,delta01]]))^(-1);
H0:=Matrix([[1,1],[delta00,delta01]]);
H1:=(Matrix([[delta21,delta31],[delta22,delta32]]))^(-1);
H2:=Matrix([[delta41,delta51],[delta42,delta52]]);
H3:=(Matrix([[delta61,delta71],[delta62,delta72]]))^(-1);
H4 := Matrix([[Y11, Y12], [Y21, Y22]]);
H5 :=( Matrix([[G11, G12], [G21, G22]]))^(-1);
H6 := Matrix([[W11, W12], [W21, W22]]);
H7 := (Matrix([[Z11, Z12], [Z21, Z22]]))^(-1);
H8 := Matrix([[V11], [V21]]);

Ht:=Matrix([[1],[delta11]]);

e1:=Matrix([[evalf(xi11(tau0)),evalf(-xi12(tau0))],[evalf(eval(diff(xi11(tau),tau),tau=tau0)),evalf(eval(diff(-xi12(tau),tau),tau=tau0))]]);

e2:=(Matrix([[evalf(xi11(tau0)),evalf(xi12(tau0))],[evalf(eval(diff(xi11(tau),tau),tau=tau0)),evalf(eval(diff(xi12(tau),tau),tau=tau0))]]))^(-1);

EE:=evalm(e1.e2);

ee1:=Matrix([[evalf(xi21(tau0)),evalf(xi22(tau0))],[evalf(eval(diff(xi21(tau),tau),tau=tau0)),evalf(eval(diff(xi22(tau),tau),tau=tau0))]]);

ee2:=(Matrix([[evalf(xi21(tau0)),evalf(-xi22(tau0))],[evalf(eval(diff(xi21(tau),tau),tau=tau0)),evalf(eval(diff(-xi22(tau),tau),tau=tau0))]]))^(-1);
EEE:=evalm(ee1.ee2);

MM:=evalm(Hi.EE.EEE.H0.H1.H2.H3.H4.H5):
R:=evalf(MM[2,1]/MM[1,1]):
Kr:=abs(evalf(R)):
KR[N]:=abs(Kr);

end do:


N:='N':
seq(KR[N],N=1..Num);
with(plots):
listplot([seq([0.5+2*(N2-1)/(Num-1),KR[N2]],N2=1..Num)]);
 

Please Wait...