clear d Dd c r f xold x xcopy ni ndf h hh ad al ar adcopy b bcopy
clear uold u res globalerr ok ans i uex Duex exacterror
if run==0
%close
figure('Position',[500 40 500 850])
else
clf
end
parameters
%
[ni,x]=pred(xmin,xmax,force,errnorm,tol);
%
if exactflag=='y';
xcopy=x;
for i=1:201;
x=xmin+(i-1)*(xmax-xmin)/200;
uexplot(i)=eval(uexact);
end
xplot=[xmin:(xmax-xmin)/200:xmax];
subplot(311),plot(xplot,uexplot,':')
figure(gcf)
title('exact solution')
x=xcopy;
end
%
%main loop
%
disp('entering main loop')
ok=-1;
%run=0;
xold=[xmin,xmax];
uold=[0,0];
while ok<0
%
%creating data vectors
%
xcopy=x;
h=x(2:ni+1)-x(1:ni);
for i=1:ni;
x=xcopy(i)+h(i)/2;
d(i)=eval(diffusion);
c(i)=eval(convection);
r(i)=eval(absorbtion);
f(i)=eval(force);
x=x+0.01*h(i);
dplus=eval(diffusion);
Dd(i)=(dplus-d(i))/(0.01*h(i));
end
%
%plotting current mesh
%
x=xcopy;
hh=[h(1),h,h(ni)];
hh=(hh(2:ni+2)+hh(1:ni+1))/2;
subplot(313),plot(x,hh)
figure(gcf)
if run==0;
title('initial meshsize');
elseif run>0
title('new meshsize')
end
pause(3);
%
%computing stiffness matrix (diagonal, left and right subdiagonals) and r.h.s.
%
[ad,al,ar]=stiff(ni,x,d,Dd,c,r,k,bctype);
b=loads(ni,x,d,Dd,c,r,f,g,bctype);
%
%solving tridiagonal system of equations
%
u=0;
adcopy=ad;
bcopy=b;
ndf=ni-1+sum(bctype);
for i=2:ndf;
ad(i)=ad(i)-ar(i-1)*al(i-1)/ad(i-1);
b(i)=b(i)-b(i-1)*al(i-1)/ad(i-1);
end
u(ndf)=b(ndf)/ad(ndf);
for i=ndf-1:-1:1;
b(i)=b(i)-b(i+1)*ar(i)/ad(i+1);
u(i)=b(i)/ad(i);
end
ad=adcopy;
b=bcopy;
%
%adding given boundary values
%
if bctype(1)==0
u=[g(1),u];
end
if bctype(2)==0
u=[u,g(2)];
end
u=u'
u=u';
%
%plotting current solution
%
subplot(313),title('current meshsize')
if exactflag=='y'
subplot(311),plot(xplot,uexplot,':'),hold on,title('exact solution');
pause(1)
end
%subplot(313),title('meshsize')
subplot(311),plot(x,u,'y-')
title('current solution')
pause(2);
hold off
%
%testing if solution ok
%
[ok,globalerr,res]=test(ni,x,u,xold,uold,d,Dd,c,r,f,errnorm,tol,ad,al,ar,bctype);
%
%refining the mesh if error too large
%
if ok<0
disp('error too large')
error=globalerr
disp('refining the mesh')
%
if exactflag=='y';
exacterror=exacterr(ni,x,u,uexact,errnorm)
end
%
pause(1)
%
[ni,x]=corr(ni,x,res,errnorm,tol);
run=run+1;
%
%accepting current solution if error ok
%
elseif ok>0
disp('error ok');
error=globalerr
%
if exactflag=='y';
exacterror=exacterr(ni,x,u,uexact,errnorm)
if errnorm=='en'&exacterror<0.1*globalerr;
disp('I notice that the exact error is much smaller than the estimated one')
disp('Probably the (midpoint) quadrature rule has underestimated')
disp('the true error !')
end
if exacterror>globalerr;
disp('I notice that the exact error is larger than the estimated one')
disp('Please check that your exact solution is correctly specified')
disp('If so, please notify the author of the code')
end
end
%
pause(1)
subplot(311),plot(x,u,'g-')
title('tested and accepted solution')
%title('solution')
pause(3);
end %test loop
end %main loop





