sundialsTB > cvodes > cvbx.m
cvbx

## PURPOSE

CVBX - CVODES example problem (serial, band)

## SYNOPSIS

This is a script file.

## DESCRIPTION

```CVBX - CVODES example problem (serial, band)
The following is a simple example problem with a banded Jacobian,
with the program for its solution by CVODES.
The problem is the semi-discrete form of the advection-diffusion
equation in 2-D:
du/dt = d^2 u / dx^2 + .5 du/dx + d^2 u / dy^2
on the rectangle 0 <= x <= 2, 0 <= y <= 1, and the time
interval 0 <= t <= 1. Homogeneous Dirichlet boundary conditions
are posed, and the initial condition is
u(x,y,t=0) = x(2-x)y(1-y)exp(5xy).
The PDE is discretized on a uniform MX+2 by MY+2 grid with
central differencing, and with boundary values eliminated,
leaving an ODE system of size NEQ = MX*MY.
This program solves the problem with the BDF method, Newton
iteration with the CVBAND band linear solver, and a user-supplied
Jacobian routine.
It uses scalar relative and absolute tolerances.
Output is printed at t = .1, .2, ..., 1.
Run statistics (optional outputs) are printed at the end.

## CROSS-REFERENCE INFORMATION

This function calls:
• CVode CVode integrates the ODE.
• CVodeFree CVodeFree deallocates memory for the CVODES solver.
• CVodeGetStats CVodeGetStats returns run statistics for the CVODES solver.
• CVodeMalloc CVodeMalloc allocates and initializes memory for CVODES.
• CVodeMonitor CVodeMonitor is the default CVODES monitoring function.
• CVodeSetOptions CVodeSetOptions creates an options structure for CVODES.
• cvbx_J CVBX_J - Jacobian functino for the CVBX example problem.
• cvbx_f CVBX_F - RHS function for the CVBX example problem
• cvbx_q CVBX_Q - quadrature function for the CVBX example problem.
This function is called by:

## SOURCE CODE

```0001 %CVBX - CVODES example problem (serial, band)
0002 %   The following is a simple example problem with a banded Jacobian,
0003 %   with the program for its solution by CVODES.
0004 %   The problem is the semi-discrete form of the advection-diffusion
0005 %   equation in 2-D:
0006 %      du/dt = d^2 u / dx^2 + .5 du/dx + d^2 u / dy^2
0007 %   on the rectangle 0 <= x <= 2, 0 <= y <= 1, and the time
0008 %   interval 0 <= t <= 1. Homogeneous Dirichlet boundary conditions
0009 %   are posed, and the initial condition is
0010 %      u(x,y,t=0) = x(2-x)y(1-y)exp(5xy).
0011 %   The PDE is discretized on a uniform MX+2 by MY+2 grid with
0012 %   central differencing, and with boundary values eliminated,
0013 %   leaving an ODE system of size NEQ = MX*MY.
0014 %   This program solves the problem with the BDF method, Newton
0015 %   iteration with the CVBAND band linear solver, and a user-supplied
0016 %   Jacobian routine.
0017 %   It uses scalar relative and absolute tolerances.
0018 %   Output is printed at t = .1, .2, ..., 1.
0019 %   Run statistics (optional outputs) are printed at the end.
0020 %
0022
0024 % Copyright (c) 2005, The Regents of the University of California.
0025 % \$Revision: 1.3 \$Date: 2006/03/07 01:19:54 \$
0026
0027 xmax = 2.0;
0028 ymax = 1.0;
0029 mx = 10;
0030 my = 5;
0031
0032 rtol = 0.0;
0033 atol = 1.0e-5;
0034
0035 t0 = 0.0;
0036 dtout = 0.1;
0037 nout = 10;
0038
0039 dx = xmax/(mx+1);
0040 dy =  ymax/(my+1);
0041
0042 % Problem data structure
0043 data.xmax = xmax;
0044 data.ymax = ymax;
0045 data.mx = mx;
0046 data.my = my;
0047 data.dx = dx;
0048 data.dy = dy;
0049 data.hdcoef = 1.0/dx^2;
0050 data.hacoef = 0.5/(2.0*dx);
0051 data.vdcoef = 1.0/dy^2;
0052
0053 % initial conditions for states
0054 t = t0;
0055 u = zeros(mx*my,1);
0056 for j = 1:my
0057   y = j * data.dy;
0058   for i = 1:mx
0059     x = i * data.dx;
0060     u(j + (i-1)*my) = x*(xmax-x)*y*(ymax-y)*exp(5.0*x*y);
0061   end
0062 end
0063
0064 % Initial condition for quadrature variable
0065 q = 0.0;
0066
0067 options = CVodeSetOptions('RelTol',rtol, 'AbsTol',atol);
0068 options = CVodeSetOptions(options,...
0069                           'LinearSolver','Band',...
0070                           'JacobianFn',@cvbx_J,...
0071                           'UpperBwidth',my,...
0072                           'LowerBwidth',my);
0073 options = CVodeSetOptions(options,...
0077
0078
0079
0080 mondata.grph = false;
0081 options = CVodeSetOptions(options,...
0082                           'MonitorFn',@CVodeMonitor,...
0083                           'MonitorData',mondata);
0084
0085 CVodeMalloc(@cvbx_f, t, u, options, data);
0086
0087
0088 ff=figure;
0089 hold on;
0090 box
0091
0092 umax = norm(u,'inf');
0093 uavg = cvbx_q(t,u,data);
0094 fprintf('At t = %f   max.norm(u) = %e\n',t, umax);
0095
0096 for i = 1:nout
0097
0098   t_old = t;
0099   uavg_old = uavg;
0100
0101   tout = t + dtout;
0102   [status,t,u, q] = CVode(tout, 'Normal');
0103
0104   if status ~= 0
0105     return
0106   end
0107
0108   uavg = cvbx_q(t,u,data);
0109   umax = norm(u,'inf');
0110   fprintf('At t = %f   max.norm(u) = %e\n',tout, umax);
0111
0112   figure(ff);
0113   plot([t_old t],[uavg_old uavg]);
0114   plot([t0 tout], [q q]/(tout-t0), 'r-');
0115   plot([tout tout], [0 q]/(tout-t0), 'r-');
0116
0117 end
0118
0119 si= CVodeGetStats
0120
0121 CVodeFree;```
Generated on Thu 05-Oct-2006 15:56:58 by m2html © 2003