We solve the following PDE in polar coordinates \(u=u(r,q)\):
\begin{align} (r\,u_{r})_{r} +\frac{1}{4r}\,(u_{q})_{q}=0,\quad 1<r<3,\quad 0<q<\pi \end{align}
with the boundary conditions \(u(1,q)=\sin(2\,q)\), \(u(3,q)=0\), and \(u(r,0)=u(r,\pi)=0\).
Thanks to Rene at reneif07@gmail.com
 |  |
download example: 2DLPSemiAnnulus.edp or return to 2D examples
//We solve the following PDE in polar coordinates
// (r*u_r)_r +(0.25/r)*(u_q)_q=0, for 1<r<3, 0<q<pi
//Boundary conditions u(1,q)=sin(2*q), u(3,q)=0, u(r,0)=u(r,pi)=0
load "Element_P3";
real beta=pi/2, r1=1, r2=3,a=1./4.;
border L1(t=-beta,beta){ x=r1*sin(t); y=r1*cos(t);label=4;}
border L2(t=-beta,beta){ x=-r2*sin(t); y=r2*cos(t);label=3;}
border L3(t=0,1){x=-r2*sin(beta)+(-r1*sin(beta)+r2*sin(beta))*t;
y=r2*cos(beta)+(r1*cos(beta)-r2*cos(beta))*t;label=2;}
border L4(t=0,1){x=r1*sin(beta)+(r2*sin(beta)-r1*sin(beta))*t;
y=r1*cos(beta)+(r2*cos(beta)-r1*cos(beta))*t;label=1;}
for (int i=1; i<=40; i++) {
int n=5*i;
mesh T0h = buildmesh(L1(n)+L2(n)+L3(n)+L4(n));
//plot(T0h,wait=1);
fespace V0h(T0h,P3);
V0h u,v,r,Ur,Uq,q;
u=0;v=0;r=sqrt(x^2+y^2);q=atan2(y,x);
Uq=sin(2.*q); Ur=((9-(r)^(2.))/(8.*r));
//
macro Grad(u) [(dx(u)*cos(q)+dy(u)*sin(q))/sqrt(r),
(dy(u)*cos(q)-dx(u)*sin(q))*sqrt(a)/sqrt(r)] // EOM
macro Lap(u,v) (Grad(u)'*Grad(v)) //
macro Jac()(r) //Jacobian
//
problem static(u,v)=int2d(T0h)(Lap(u,v)*Jac)
+ on(1,2,3,u=0)
+ on(4,u=Uq);
//exact solution
V0h f=Uq*Ur;
//
static;
V0h err=u-f;
V0h Int=(err)^2.*Jac;
real l2e=sqrt(int2d(T0h)(Int));
cout << " L2-error " << " = "<<l2e<< endl;
plot(u,wait=1,value=1,fill=1,cmm="Max(u_numeric)="+u[].max+", Min(u_numeric)="+u[].min+", L^2-error="+l2e);
plot(f,value=1,wait=1,fill=1,cmm="Max(u_exact)="+f[].max+", Min(u_exact)="+f[].min+", L^2-error="+l2e);
}