# sandbox/easystab/differential_equation.m

# Solving a linear differential equation

This code is built upon diffmat_dif1D.m where the differentiation matrices are built and validated. Here we solve a linear non-homogeneous differential equation in 1D with non-homogeneous boundary conditions.

We show here how to use the differentiation matrix to code the differential equation as a linear system and impose Dirichlet and Neuman boundary conditions.

To learn more about differentiation matrices, see pedagogy#differentiation-matrices and for boundary conditions, please see pedagogy#boundary-conditions.

```
clear all; clf
% parameters
L=2*pi; % domain length
N=25; % number of points
% 1D differentiation matrices
[D,DD,wx,x]=dif1D('fd',0,L,N,3);
```

# The differential equation

We solve this \displaystyle f_{xx}=1 with the boundary conditions \displaystyle f(1)=0, f_x(L)=1

So you see that this is a non-homogeneous equation, since there is a forcing term. And we have both Dirichlet condition (at x=0) and Neumann boundary condition (at x=L). You see that we have used the first order differentiation matrix also when imposing the boundary condition. Off course this is very natural since these matrices are the way we have transformed the operation of differentiation for our discrete systems.

What we have done is the following. We have replaced the first equation (the first line of the system matrix *A* by a new equation that tells that the value of *f* at the first gridpoint should be equal to 0. We did this using the first line of the identity matrix *I*. Indeed, this first line multiplied by the vector of *f* is the first value in *f*. We then did the same thing for the Neumann boundary condition at the last gridpoint. We have replaced the line in *A* by the last line of the differentiation matrix *D*. Indeed, this last line multiplied by the vector of *f* will give you the derivative at L.

```
% I build the linear differential equation
A=DD;
% boundary conditions
I=eye(N);
A([1,N],:)=[I(1,:); D(N,:)];
b=1+zeros(N,1); b([1,N])=[0,1];
% solve the system
f=A\b;
% plotting
plot(x,f,'b.-',x,x.^2/2+(1-L)*x,'r-');
xlabel('x');ylabel('f');
legend('numerical','theory')
xlim([0,L]); grid on
set(gcf,'paperpositionmode','auto')
print('-dsvg','-r75','plot.svg')
```

# Links

We do more or less the same things but in 2D in poisson2D.m and in 3D in poisson3D.m.

In variable_coef.m we solve in 1D a differential equation with variable coefficients.

# Exercices/Contributions

- Please put the two boundary conditions on the same side of the domain (for instance, impose the value of f and its derivative at x=0) –> differential_equation_sameside.m
- Please solve f_{xx}=\sin(x) –> differential_equation_sin.m
- Please check the solution and convergence of the system a
*Df+b*DDf=c, where a,b stands for scalars and c a vector –> differential_equation_secondorder.m - Please check the convergence with the number of gridpoints –> convergence_of_the_solution.m
- Please use Chebychev differentiation matrix chebdif.m and compare with the finite differences -> differential_equation_chebychev.m
- Please use a non-constant coefficient i your differential equation a(x)f_{xx}=0 (first find the way to code that using differentiation matrices) *Please check convergence of the system aDDf+bDf+cf+d=0 classical_differential_equation.m