I am trying to solve a pair of coupled linear PDEs using backward Euler and Newton-Raphson

I am tring to solve a coupled pair of 1D PDEs in time and space using the backward Euler and Newton-Raphson. One of them is a hyperbolic equation and the other one is a parabolic equation. The method I am using is this:

1) I write the pair of variables at a given time in one column vector y=[rho; v;]

2) The backward Euler method gives me a set of linear equations in the components of y.

3) I use the Newton-Raphson method to solve the set of linear equations.

4) I use the solution for the previous timestep as an initial guess for the current timetep.

I tried this method for a simple nonlinear ODE, y'=1-y^2, and it worked beautifully. The scheme worked very quickly, only needing 1 interation in the NR scheme. When I try it on my linear scheme, it blows up, and I am at a loss really.

The equations I am solving are:

r_t+0.1*v_x=0

v_t=0.01*v_xx

The code I write is:

%This is to test the method of solution of a ODE via Newton-Raphson method.
clear;clc;
N=200;
M=1000;
t=linspace(0,5,M);
x=linspace(0,1,N); %x-axis variables
dx=x(2);
dt=t(2);
D=0.01;
alpha=dt/dx;
beta=D*dt/dx^2;
y=zeros(2*N,1); %Solution variables;
tol=1e-8;
%Add in initial conditions
y_old(1:N)=1;
y_old(N+1:2*N)=0;
rho=zeros(M,N);
v=zeros(M,N);
rho(1,:)=y_old(1:N);
y_old=[rho(1,:)'; v(1,:)'];
%Compute the solution
for i=2:N
    y_test=1.001*y_old;
    y=Newton(tol,alpha,beta,dx,t(i),y_old,y_test);
    y_old=y;
    rho(i,:)=y(1:N)';
    v(i,:)=y(N+1:2*N)';
end

function [J] = jacobian(alpha,beta,dx,t,y_old,y)
% Computes the Jacobian matrix of the function f(x)
% f: function handle that returns a vector of function values
% x: column vector of independent variables

% Number of function values and independent variables
n = length(y); % number of independent variables

% Initialize Jacobian matrix
eps=1e-8; % could be made better
J = zeros(n,n);
T=zeros(n,1);
for j=1:n
    T(j)=1;
    fPlus = test_eqns(alpha,beta,dx,t,y_old,y+eps*T);
    fMinus = test_eqns(alpha,beta,dx,t,y_old,y-eps*T);
    J(:,j) = (fPlus-fMinus)/(2*eps);
    T(j)=0;
end
end

% Define the Modified Newton-Raphson method
function [y] = Newton(tol,alpha,beta,dx,t,y_old,y0)
    iter = 0;
    maxiter = 10;
    error=10;
    while (error > tol) && (iter < maxiter)
        J0 = jacobian(alpha,beta,dx,t,y_old,y0);
        f=-test_eqns(alpha,beta,dx,t,y_old,y0)';
        delta_y = (f/J0)';
        y_new = y0 + delta_y;
        y0 = y_new;
        error=norm(f);
        iter = iter + 1;
    end
    y=y0;
end

function f=test_eqns(alpha,beta,dx,t,y_old,y0)

N=length(y0);
n=floor(N/2);
f=zeros(N,1);
for i=2:n-1
    %Conservation of mass
    f(i)=y0(i)-y_old(i)-alpha*1.5*(y0(n+i)-y0(n+i-1));
    %Conservation of momentum
    f(n+i)=beta*y0(i+1)-(1+2*beta)*y0(i)+beta*y0(i-1)+y_old(i);
end

%Add in the boundary conditions
%Conservation of mass
f(1)=y0(2)-y0(1);
f(n)=y0(n)-y0(n-1);
%Conservation of momentum.
f(n+1)=y0(n+1);
f(N)=y0(N)-y0(N-1)-dx*t.^2.*exp(-0.1*t);
end

I really have no idea what is going on. I have used an upwind scheme for the hyperbolic equation, as that was advised from people in the know.

[Torsten]

As the second equation is an easy-to-solve parabolic PDE and independent of the first, I suggest you try first to solve it independently of the first equation. This should be a good test for your code from above although I'd suggest to use the method-of-lines in combination with ODE15S to solve the resulting system of ordinary differential equations. But if you want, you can of course stick to your code above and program the time integration as Implicit Euler on your own.

The second equation looks strange to me: If r(t,x) is a solution, so is r(t,x) + F(x) for an arbitrary function F(x) since d/dt (F(x)) = 0. So you will have to prescribe r(0,x). But this is the only condition you can impose since the complete behaviour over t follows from v. Setting

%Conservation of mass
f(1)=y0(2)-y0(1);
f(n)=y0(n)-y0(n-1);

is thus not possible for rho.

[Matthew Hunt]

They're connected. It's a prototype sytem to the one I really want to solve,so ad hoc methods are not the sort of thing that I would like to consider.

I'm unwilling to use ODE15, as it hasn't worked in the past, and I would rather not waste any more time trying to get it to work.

The second equation is a diffusion equation.

I don't understand what you're attempting to say.

[Torsten]

They're connected. It's a prototype sytem to the one I really want to solve,so ad hoc methods are not the sort of thing that I would like to consider.

What do you mean by ad hoc methods in this context ?

I'm unwilling to use ODE15, as it hasn't worked in the past, and I would rather not waste any more time trying to get it to work.

ODE15s is the well-tested alternative to your Implicit Euler + Newton's method. So if your discretized equations work with Implicit Euler + Newton's method, they should equally work with ODE15S.

I don't understand what you're attempting to say.

I attempt to say that you should test your code with the mathematically simple diffusion equation first before you try to solve both equations simultaneously. I see a coupling of the first equation to the second, but not a coupling of the second to the first - at least in the section where you wrote the equations in a mathematical notation.

[Matthew Hunt]

Ad hoc methods being methods that only work with one particular system of equations and on in a general case.

I could never get ODE5 to work, and it's a black box where you can't really see what's going on. The method that I'm working on, you can see all the gory details.

My code works with the viscous berguers equation just fine.

Answer 1

[Garmit Pant]

Hello Matthew

I understand that you are trying to solve a coupled pair of PDEs using a combination of backward Euler and modified Newton-Raphson method.

MATLAB does not have built-in functions for backward Euler and Newton-Raphson method. You can consult the following links that have implementations of these methods:

1. Multiple implementations of Backward Euler method that solves the backward Euler equation are available on the internet. You can consult any of the sources available online to check the implementation.
2. https://in.mathworks.com/help/optim/ug/fsolve.html - The MATLAB function ‘fsolve’ can be used to solve a system of non-linear equations (Requires Optimazation Toolbox)
3. https://www.mathworks.com/matlabcentral/fileexchange/68885-the-newton-raphson-method - Implementation of the ‘Newton - Raphson Method’, to solve various polynomial and transcendental equations

You can also try to modify the boundary conditions and set the initial conditions accordingly.

MATLAB offers the function ‘pdepe’ to solve 1-D parabolic and elliptic PDEs. The function also produces results for hyperbolic equations. You can follow the example given in the following documentation page and try solving the equations using the function ‘pdepe’:

https://www.mathworks.com/help/matlab/ref/pdepe.html

I hope this helps!

Best Regards

Garmit