Code has numerical issues that I don't understand

I managed to get a system of PDEs down to a large system of ODEs via the method of lines. I managed to code it up using ode45 but I still have numerical errors when the code is run, and I don't quite understand why. The basic physics of the model is that of a porous 1D bar where the porosity and length decreases, whereas the velocity and temperature increases. I think I have reasonable initial conditions and I've checked my discritisation and it looks okay. I get te error:

Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN.

> In odemassexplicit>ExplicitSolverHandleMass34 (line 67)

In odefcncleanup>@(t,y)oldFcn(t,y,inputArgs{:}) (line 11)

In ode45 (line 315)

In sintering (line 12)

clear;clc;
N=10;
c=1;
h=1/N;
tspan=[0 20];
y_0=zeros(1,3*N+1);
y_0(1:N)=1;
y_0(N+1:2*N)=0;
y_0(2*N+1:3*N)=1;
y_0(3*N+1)=1;
opts = odeset('Mass',@(t,y) mass(y,h,c));
[t,y] = ode45(@(t,y) sint(N,y),tspan,y_0,opts);
function dydt = sint(N,y)
%This is the RHS of the sintering equations
g=9.81;
alpha=0.1;
h=1/N;
kappa=1;
x=linspace(0,1,N);
dydt=zeros(3*N+1,1);
dydt(3*N+1)=y(2*N);
%These are the bulk equations
for i=2:N-1
    dydt(i)=y(2*N)*(y(i+1)-y(i-1))/(2*h*y(3*N+1))-(y(3*N+1)/(2*h))*(y(i+1)*y(N+1+1)-y(i-1)*y(N+i-1));
    A=(y(2*N)*y(i)/(2*h*y(3*N+1)))*(y(i+1)-y(i-1))-(y(i)*y(N+i)/(2*h*y(3*N+1)))*(y(N+i+1)-y(N+i-1))+1/(2*h*y(3*N+1))*(P_L(y(i+1))-P_L(y(i-1)));
    C=(1/(h^2*y(3*N+1)^2))*0.5*(zeta(y(i+1))+zeta(y(i)))*y(N+i+1)-0.5*(zeta(y(i+1))+zeta(y(i))+zeta(y(i))+zeta(y(i-1)))*y(N+i)+0.5*(zeta(y(i+1))+zeta(y(i-1)))*y(N+i-1);
    D=-g+alpha/(2*h*y(3*N+1))*(y(2*N+i+1)-y(2*N+i-1));
    dydt(N+i)=A+C+D;
    E=-((y(N+i)/y(3*N+1))-y(2*N)/y(3*N+1).^2)*((y(N+i+1)-y(N+i-1))/(2*h));
    F=(kappa*y(i)/(h^2*y(3*N+1)))*(y(2*N+i+1)-2*y(2*N+i)+y(2*N+i-1))-(P_L(y(i))/y(3*N+1))*((y(N+i+1)-y(N+i-1))/(2*h));
    G=(alpha*y(i)*y(2*N+i)/h^2)*(y(N+i)/y(3*N+1)^2-y(2*N)/y(3*N+1))*(y(N+i+1)-2*y(N+i)+y(N+i-1));
    dydt(2*N+i)=E+F+G;
end

%This is the boundary condition for the density;
dydt(1)=0;
dydt(N)=y(2*N)*(y(N)-y(N-1))/(h*y(3*N+1));

%This is the boundary condition for the velocity

dydt(N+1)=(P_L(y(2))-P_L(y(1)))/(h*y(3*N+1))+(y(N+2)*(zeta(y(2))-zeta(y(1))))/(h^2*y(3*N+1)^2)-g+alpha*(y(2*N+2)-y(2*N+1))/(h*y(3*N+1));
A=y(2*N)*y(N)*(y(N)-y(N-1))/(h*y(3*N+1))+2*y(N)*y(2*N)*(P_L(y(N))+alpha*T_a)/(y(3*N+1)*zeta(y(N)))+(P_L(y(N))-P_L(y(N-1)))/(h*y(3*N+1));
B=(0.5*(3*zeta(y(N))-zeta(y(N-1)))*(y(2*N-1)-2*h*(P_L(y(N)))+alpha*T_a)/zeta(y(N))+0.25*y(2*N)*(7*zeta(y(N))-zeta(y(N-1)))+0.5*y(2*N-1)*(zeta(y(N))+zeta(y(N-1))))/(h^2*y(3*N+1)^2);
C=-g+alpha*(y(3*N)-y(3*N-1))/(h*y(3*N+1));
dydt(2*N)=A+B+C;

%This is the boundary condition for Temperature

dydt(2*N+1)=y(1)*(kappa*y(3*N+1)-2*kappa*T_a+kappa*(2*T_a-y(3*N+1)))/(h^(2)*y(3*N+1)^2)-P_L(y(1))*y(N+2)/h;
dydt(3*N)=y(N)*(kappa*(2*T_a)-y(3*N-1)-2*kappa*T_a+kappa*y(3*N-1))-P_L(y(N))*(y(2*N)-y(2*N-1))/(h*y(3*N+1));
end

function M=mass(y,h,c)
N=length(y);
n=floor((N-1)/3);
l=ones(1,n);
M=zeros(N,N);
%Insert the conservation of mass terms
M(1:n,1:n)=diag(l,0);

%Insert the coefficients for the conservation of momentum
M(n+1:2*n,n+1:2*n)=diag(y(1:n),0);

%Insert the coefficients for the conservation of energy
M(2*n+1:3*n,2*n+1:3*n)=diag(c*y(1:n),0);  %Diagonal elements
%Compute the off diagonal elements
l_sub=ones(1,n-1)/(2*h);
M_sub=diag(l_sub,-1)-diag(l_sub,1);
M_sub(1,1)=1/h;
M_sub(1,2)=-1/h;
M_sub(n,n-1)=1/h;
M_sub(n,n)=-1/h;
M((2*n+1):3*n,(n+1):2*n)=M_sub;
M(N,N)=1;
end

function y=P_L(x)
gamma=1;
r_0=1;
y=3*gamma*(1-x).^2/r_0;
end

function y=zeta(x)
eta_0=1;
y=2*(1-x).^3*eta_0/(3*x);
end

function y=T_a(t)
y=1;
end

Can you suggest might be wrong?

[Torsten]

How should we be able to suggest something if we do not even know the PDE you tried to discretize ?

[Matthew Hunt]

I don't think that's the issue. The system of PDEs was discretised using finite differences. That's good enough for most equations.

[Bjorn Gustavsson]

@Matthew Hunt: it might not be the issue you have problems with (but also remember that the phrase "famous last words..." exists for a reason). But it is the issue for other peoples possibilities to give advice. You need to make it easy for others to understand what the code does - fewer will have the enthusiasm to read your code and try to deduce what it is suppoesed to solve...

[Matthew Hunt]

I had a system of three PDEs in density, velocity and temperature. I then discritised in space via finite differences to obtain a system of ODEs. This allows me to use ode45 to solve the lot. I couldn't use the PDE solver because of the mixed derivative.

[Torsten]

I had a system of three PDEs in density, velocity and temperature. I then discritised in space via finite differences to obtain a system of ODEs. This allows me to use ode45 to solve the lot. I couldn't use the PDE solver because of the mixed derivative.

We know this, but obviously your discretization causes trouble. And discretizing mixed derivatives together with boundary conditions is not an easy task.

So I think it's natural to ask whether your discretization is adequate for the terms in your PDEs you are trying to solve. And in order to answer this, we have to see your PDEs together with initial and boundary conditions. But - as said - the problem seems to be hard and there is no guarantee we can help. So if it's too much effort for you to supply the problem, we are not offended :-)

[Bjorn Gustavsson]

So coupled fluid equations? Continuity equation, momentum equation and energy equation?

Like

and so on? Please give us those equations to start looking at...

[Torsten]

@Bjorn Gustavsson

It seems to be a continuation of this question:

https://de.mathworks.com/matlabcentral/answers/1891380-solving-pdes-with-mixed-derivatives?s_tid=srchtitle

I don't know where the mixed derivatives in the energy equation stem from.

[Matthew Hunt]

I derived the temperature equation via thermodynamics and as a result obtained a term with a mixed derivative. I've not seen it before and thought it was an interesting inclusion into the model. Physically I'm tring to model the temperature equation for sintering, my general approach is taken from this paper:

http://gitstech.weebly.com/uploads/1/0/0/2/10023880/theory_of_sintering_.pdf

I have 3 basic PDEs:

Conservation of mass:

Newton's second law:

And the temperature equation:

I discritised these equations in space and put them into ode45.

[Torsten]

In the momentum equation, I don't understand d(nu)/du and dD/d(epsilon).

In the energy equation, I don't understand d^2(nu)/(dx dx).

Apart from t and x, are u and epsilon further independent variables in the direction of which you have to discretize ?

Where is the mixed derivative ?

If you really want to solve continuity and momentum equation, you have to apply very specialized discretization schemes that were developed for fluid dynamics applications. You will find code for them under

https://www.clawpack.org/

Under the examples section, look up "Euler equations".

[Matthew Hunt]

The d/du was a typo, and was corrected as was the mixed derivative. The mixed derivative comes from including the thermal stress tensor within the derivation of the temperature equtation. Here epsilon is the rate of strain tensor in 1D, and D is the dissipation potential and can be considered as a function of \partial_{x}v,, \rho and \partial_{x}\rho, it's equal to:

and zeta is the bulk moluli.

I don't want to use clawpack. I know how to discretise these equations, they should work.

[Torsten]

I know how to discretise these equations, they should work.

Believe me: these equations are very difficult to solve and adequate discretization schemes fill books in the literature.

One is tempted to use standard discretization schemes, but they won't work.

As long as velocity and pressure are given, standard discretization schemes work (e.g. for species or temperature). As soon as you try to add the continuity and momentum equations with a naive approach, the results become unphysical.

But the best way to learn is to fail - and finally to know why one failed.