|
>
|
restart:with(plots):with(PDETools):with(LinearAlgebra):
m:=(100.+(t/To)^(1/10))*(t/To)^(2/5)/(100.+(t/To)^(1/2)); To:=1;
lame_cyl:=r^2*diff(u(r,t),r$2)+r*diff(u(r,t),r)+m*u(r,t)=0; bound_con:=u(1,t)=1,u(2,t)=-1; a:=1; b:=2;
N:=20; h:=(b-a)/N;R:=k->a+k*h; U1:=(k,t)->(u[k+1,t]-u[k-1,t])/(2*h); U2:=(k,t)->(u[k+1,t]-2*u[k,t]+u[k-1,t])/h^2;
e[0]:=u[0,t]=rhs(bound_con[1]);e[N]:=u[N,t]= rhs(bound_con[2]);
for k from 1 to N-1 do e[k]:=eval(lame_cyl,{r=R(k),u(r,t)=u[k,t], diff(u(r,t),r)=U1(k,t),diff(u(r,t),r$2)=U2(k,t)} );od;
row[0]:=[rhs(bound_con[1]),seq(0,j=1..N-1)];
row[1]:=[coeff(lhs(e[1]),u[0,t]),coeff(lhs(e[1]),u[1,t]),coeff(lhs(e[1]),u[2,t]),seq(0,j=1..N-3)];
for n from 2 to N-1 do row[n]:=[seq(0,j=1..n-2),coeff(lhs(e[n]),u[n-1,t]),coeff(lhs(e[n]),u[n,t]),coeff(lhs(e[n]),u[n+1,t]),seq(0,j=1..N-1-n)];od;
row[N]:=[seq(0,j=1..N-1),rhs(bound_con[2])];
row_matrix:=[row[0],seq(row[n],n=1..N-2),row[N]];A:=(row_matrix);
U:=Vector([seq(u[j,t],j=1..N)]);B:=Vector([rhs(bound_con[1]),seq(rhs(e[j]),j=1..N-2),rhs(bound_con[2])]);
with(linalg):A.U = B;;;A1:=inverse(A);U:=linsolve(A,B);
result:=[seq(evalf(subs(t=10,U(k,t)),1),k=0..N)]:
print("Ham epsilon[theta](1,t)");plot3d(-U[N-1]/r,r=1..1.000001,t=0..100
);print("Ham epsilon[theta](2,t)");plot3d(U[N-1]/r,r=2..2.000001,t=0..100);
u(r,t):=U[N-1];;Ee:=0.5;;epsilon[theta](r,t):=u(r,t)/r;E[theta](r,t) :=
(100/(t/To)^.1+1)*Ee; sigma[theta](r,t):=epsilon[theta](r,t)*E[theta](r,t) ;Ee:=0.5;sigma[a](t):=normal(subs(r=1,-sigma[theta](r,t)));sigma[b](t):=normal(subs(r=2,sigma[theta](r,t)));;;with(plottools):with(plots):plot([sigma[a](t),sigma[a](t)],t=0..10^10,style=[point,line],symbol=cross,color=[blue,black],thickness=[1],legend=[`sigma[an](t)`,`sigma[a](t)`],axes=box,symbolsize=[10]);plot([sigma[b](t),sigma[b](t)],t=0..10^10,style=[point,line],symbol=diamond,color=[red,black],thickness=[1],legend=[`sigma[bn](t)`,`sigma[b](t)`],axes=box,symbolsize=[10]);print("HAM
u(r,t)");
smartplot3d(u(r,t));
|
