FEATool Multiphysics  v1.16.5 Finite Element Analysis Toolbox
ex_nonlinear_pde1.m File Reference

## Description

EX_NONLINEAR_PDE1 1D nonlinear PDE example.

[ FEA, OUT ] = EX_NONLINEAR_PDE1( VARARGIN ) Solves the equation (x*u*u_x)_x = 0 with Dirichlet boundary conditions u(x=x1)=u1 and u(x=x2)=u2. Compares against the analytical solution u(x) = sqrt(2*c1*log(x)+2*c2), with c1 = -(u1^2-u2^2)/ (log(x2)-log(x1))/2 and c2 = (log(x2)*u1^2-u2^2*log(x1))/(log(x2)-log(x1))/2. Accepts the following property/value pairs.

Input       Value/{Default}        Description
-----------------------------------------------------------------------------------
hmax        scalar {1/10}          Grid cell size
sfun        string {sflag1}        Finite element shape function
u1          scalar                 Dirichlet boundary value at x = x1
u2          scalar                 Dirichlet boundary value at x = x2
x1          scalar                 Left most domain x-coordinate
x2          scalar                 Right most domain x-coordinate
ischeme     scalar {-1}            Solver scheme (<0 = stationary)
iplot       scalar 0/{1}           Plot solution (=1)
.
Output      Value/(Size)           Description
-----------------------------------------------------------------------------------
fea         struct                 Problem definition struct
out         struct                 Output struct


# Code listing

 cOptDef = { ...
'hmax',     1/10;
'sfun',     'sflag1';
'u1',       1;
'u2',       2;
'x1',       3;
'x2',       5;
'ischeme', -1;
'iplot',    1;
'tol',      5e-3;
'fid',      1 };
[got,opt] = parseopt(cOptDef,varargin{:});
fid       = opt.fid;

% Grid generation.
nx       = round( 1/opt.hmax );
fea.grid = linegrid( nx, opt.x1, opt.x2 );

% Coefficients and expressions.
fea.sdim  = { 'x' };
fea.expr  = { 'u1'     opt.u1;
'u2'     opt.u2;
'x1'     opt.x1;
'x2'     opt.x2;
'c1'     '-(u1^2-u2^2)/(log(x2)-log(x1))/2';
'c2'     '(log(x2)*u1^2-u2^2*log(x1))/(log(x2)-log(x1))/2';
'u_ref'  'sqrt(2*c1*log(x)+2*c2)' };

% Physics mode and equation settings.
fea = addphys( fea, @customeqn, { 'u' } );
fea.phys.ce.sfun = { opt.sfun };
fea.phys.ce.eqn.seqn = 'u'' - x*u*ux_x = 0';
fea.phys.ce.eqn.coef = { 'u0_ce', 'u_0', 'Initial condition for u', { 'u1+(u2-u1)*(x-x1)/(x2-x1)' } };

% Boundary settings.
fea.phys.ce.bdr.coef = { 'bcnd_ce', 'Dirichlet/Neumann boundary conditions',...
'Dirichlet/Neumann boundary conditions', ...
{'Dirichlet, r_u', 'Neumann, g_u' }, ...
{ 1 1 }, [], { 'u1' 'u2' } };

% Parse and solve problem.
fea = parsephys( fea );
fea = parseprob( fea );
if( opt.ischeme<1 )
fea.sol.u = solvestat( fea, 'fid', fid, 'init', {'u0_ce'} );
else
fea.sol.u = solvetime( fea, 'fid', fid, 'init', {'u0_ce'}, 'ischeme', opt.ischeme );
end

% Postprocessing.
if( opt.iplot>0 )
figure
subplot(3,1,1)
postplot( fea, 'surfexpr', 'u', 'axequal', 'off' )
title('Solution u')
subplot(3,1,2)
postplot( fea, 'surfexpr', 'u_ref', 'axequal', 'off' )
title('Exact solution')
subplot(3,1,3)
postplot( fea, 'surfexpr', 'abs(u_ref-u)', 'axequal', 'off', 'evalstyle', 'exact' )
title('Error')
end

% Error checking.
xi  = [1/2; 1/2];
err = evalexpr0( 'abs(u_ref-u)', xi, 1, 1:size(fea.grid.c,2), [], fea );
ref = evalexpr0( 'u',            xi, 1, 1:size(fea.grid.c,2), [], fea );
err = sqrt(sum(err.^2)/sum(ref.^2));

if( ~isempty(fid) )
fprintf(fid,'\nL2 Error: %f\n',err)
fprintf(fid,'\n\n')
end

out.err  = err;
out.pass = out.err<opt.tol;
if( nargout==0 )
clear fea out
end