[hunt_mat]

In order to build up solving a coupled set of equations, I was going to solve the viscous Burger's equation:

u_t+uu_x=Dy_xx

along with the boundary conditions u(t,0)=0, and u_x(t,1)=0. I have had some success with the backward Euler method as being unconditionally stable, and I understand that I need to upwind. I wrote te following code:

%This is to test the method of solution of a ODE via Newton-Raphson method.
clear;clc;
M=3000;
N=200;
D=0.1;
x=linspace(0,1,N); %x-axis variables
t=linspace(0,20,M);
dx=abs(x(2)-x(1));
dt=t(2);
a=dt/dx;
b=D*dt/(dx^2);
u=zeros(M,N); %Solution variables;
u_0=exp(-25*(x-0.5).^2);
u(1,:)=u_0;
y_old=u_0;
tol=1e-6;
%Compute the solution
for i=2:M
    y_old=u(i-1,:);
    y_test=u(i-1,:)';
    r=Newton(tol,a,b,y_old,y_test);
    u(i,:)=r';
end

for i=1:2
    plot(x,u(i,:));
    pause(0.1);
end

function [J] = jacobian(a,b,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-5; % could be made better
J = zeros(n,n);
T=zeros(n,1);
for i=1:n
    T(i)=1;
    fPlus = vb(a,b,y_old,y+eps*T);
    fMinus = vb(a,b,y_old,y-eps*T);
    J(:,i) = (fPlus-fMinus)/(2*eps);
    T(i)=0;
end
end

% Define the Modified Newton-Raphson method
function y = Newton(tol,a,b,y_old,y0)
    iter = 0;
    maxiter = 10;
    error=10;
    while (error > tol) && (iter < maxiter)
        J0 = jacobian(a,b,y_old,y0);
        delta_y = -vb(a,b,y_old,y0)'/J0;
        y_new = y0 + delta_y';
        y0 = y_new;
        error=norm(vb(a,b,y_old, y0));
        iter = iter + 1;
    end
    y=y0;
end

function z=vb(a,b,y_old,y)
N=length(y);
z=zeros(N,1);
for i=2:N-1
    z(i)=y(i)-y_old(i)+a*y(i)*(y(i+1)-y(i))-b*(y(i+1)-2*y(i)+y(i-1));
end
%Use the boundary conditions
z(1)=y(1);
z(N)=y(N)-y_old(N)-b*(y(N-1)-y(N));
end

I can't seem to get it to work. Is there something I'm doing wrong? Are the boundary conditions incorrectly set? Are the boundary conditions incompatible?

[hunt_mat]

The issue here was the size of b was too large. I reduced the size and I got a beautiful solution which does all the funky things that it should do.

I'm curious as I thought the backward Euler was supposed to be unconditionally convergent? What's going on here?

[FMDenaro]

Originally Posted by hunt_mat

In order to build up solving a coupled set of equations, I was going to solve the viscous Burger's equation:

u_t+uu_x=Dy_xx

along with the boundary conditions u(t,0)=0, and u_x(t,1)=0. I have had some success with the backward Euler method as being unconditionally stable, and I understand that I need to upwind. I wrote te following code:

%This is to test the method of solution of a ODE via Newton-Raphson method.
clear;clc;
M=3000;
N=200;
D=0.1;
x=linspace(0,1,N); %x-axis variables
t=linspace(0,20,M);
dx=abs(x(2)-x(1));
dt=t(2);
a=dt/dx;
b=D*dt/(dx^2);
u=zeros(M,N); %Solution variables;
u_0=exp(-25*(x-0.5).^2);
u(1,:)=u_0;
y_old=u_0;
tol=1e-6;
%Compute the solution
for i=2:M
    y_old=u(i-1,:);
    y_test=u(i-1,:)';
    r=Newton(tol,a,b,y_old,y_test);
    u(i,:)=r';
end

for i=1:2
    plot(x,u(i,:));
    pause(0.1);
end

function [J] = jacobian(a,b,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-5; % could be made better
J = zeros(n,n);
T=zeros(n,1);
for i=1:n
    T(i)=1;
    fPlus = vb(a,b,y_old,y+eps*T);
    fMinus = vb(a,b,y_old,y-eps*T);
    J(:,i) = (fPlus-fMinus)/(2*eps);
    T(i)=0;
end
end

% Define the Modified Newton-Raphson method
function y = Newton(tol,a,b,y_old,y0)
    iter = 0;
    maxiter = 10;
    error=10;
    while (error > tol) && (iter < maxiter)
        J0 = jacobian(a,b,y_old,y0);
        delta_y = -vb(a,b,y_old,y0)'/J0;
        y_new = y0 + delta_y';
        y0 = y_new;
        error=norm(vb(a,b,y_old, y0));
        iter = iter + 1;
    end
    y=y0;
end

function z=vb(a,b,y_old,y)
N=length(y);
z=zeros(N,1);
for i=2:N-1
    z(i)=y(i)-y_old(i)+a*y(i)*(y(i+1)-y(i))-b*(y(i+1)-2*y(i)+y(i-1));
end
%Use the boundary conditions
z(1)=y(1);
z(N)=y(N)-y_old(N)-b*(y(N-1)-y(N));
end

I can't seem to get it to work. Is there something I'm doing wrong? Are the boundary conditions incorrectly set? Are the boundary conditions incompatible?

The velocity is always negative?


a*y(i)*(y(i+1)-y(i))

[hunt_mat]

Originally Posted by FMDenaro

The velocity is always negative?


a*y(i)*(y(i+1)-y(i))

Yes, the velocity is always negative. This is part of a model that deals with a compressible rod essentially. The compression is due to temperature and will always be towards the fixed end. the length, L, of the rod is given by dL/dt=v(t,N). For the length of the rod to decrease, the velocity must be negative.

I thought that backward Euler was unconditionally convergent?

[FMDenaro]

Originally Posted by hunt_mat

Yes, the velocity is always negative. This is part of a model that deals with a compressible rod essentially. The compression is due to temperature and will always be towards the fixed end. the length, L, of the rod is given by dL/dt=v(t,N). For the length of the rod to decrease, the velocity must be negative.

I thought that backward Euler was unconditionally convergent?

But that means you have an inflow at x=1, u(1,t) = u_inflow and an outflow at x=0 when you can use the zero derivative, right?

Why do you use Newton method?

I don’t understand the implicit part.

[hunt_mat]

Originally Posted by FMDenaro

But that means you have an inflow at x=1, u(1,t) = u_inflow and an outflow at x=0 when you can use the zero derivative, right?

Why do you use Newton method?

I don’t understand the implicit part.

I don't have any inflow/outflow for my problem. I have conservation of mass and that shrinks the domain which I fix with domain fixing. Nothing goes in or out of the region.

I use Newton-Raphson because it's the only technique I know of that I can use to solve a set of nonlinear algebraic equations.

[FMDenaro]

Originally Posted by hunt_mat

I don't have any inflow/outflow for my problem. I have conservation of mass and that shrinks the domain which I fix with domain fixing. Nothing goes in or out of the region.

I use Newton-Raphson because it's the only technique I know of that I can use to solve a set of nonlinear algebraic equations.

there is no mass in your equation, you are solving a parabolic PDE (or perturbed hyperbolic for small values of D) and you MUST prescribe the correct mathematical BCs to have a well-posed problem.
If you think that that the assumption of no-inflow and no-outflow is acceptable, use periodic conditions or set u(0,t)=u(1,t)=0. Otherwise your problem is wrong.


Now, about the implicit scheme, I suppose you solve


u(i,n+1)=u(i,n)-u(i,n+1)*(dt/dx)*(u(i+1,n+1)-u(i,n+1)) +
+D*dt/(dx^2)*(u(i+1,n+1)-2*u(i,n+1)+u(i-1,n+1))


right?


Generally, the equation is linearized, in your case you have a first order accurate scheme, therefore in the CFL number you can set u(i,n+1)=u(i,n) and solve the linear algebric system.
However, this non-conservative scheme is not the best to use.

[hunt_mat]

Originally Posted by FMDenaro

there is no mass in your equation, you are solving a parabolic PDE (or perturbed hyperbolic for small values of D) and you MUST prescribe the correct mathematical BCs to have a well-posed problem.
If you think that that the assumption of no-inflow and no-outflow is acceptable, use periodic conditions or set u(0,t)=u(1,t)=0. Otherwise your problem is wrong.


Now, about the implicit scheme, I suppose you solve


u(i,n+1)=u(i,n)-u(i,n+1)*(dt/dx)*(u(i+1,n+1)-u(i,n+1)) +
+D*dt/(dx^2)*(u(i+1,n+1)-2*u(i,n+1)+u(i-1,n+1))


right?


Generally, the equation is linearized, in your case you have a first order accurate scheme, therefore in the CFL number you can set u(i,n+1)=u(i,n) and solve the linear algebric system.
However, this non-conservative scheme is not the best to use.

There will be a mass equation, I am building up the system equation by equation. I plan to get the mass equation involved when I can get the velocity equation (given a constant density for now). The equations can't be put into conservative form I don't think.

So you're saying u(i,n+1)=u(i,n)+du, and obtain a matrix for du?

[FMDenaro]

Originally Posted by hunt_mat

There will be a mass equation, I am building up the system equation by equation. I plan to get the mass equation involved when I can get the velocity equation (given a constant density for now). The equations can't be put into conservative form I don't think.

So you're saying u(i,n+1)=u(i,n)+du, and obtain a matrix for du?

Yes, you have a low accuracy order in time, the linearization is suitable

[hunt_mat]

Originally Posted by FMDenaro

Yes, you have a low accuracy order in time, the linearization is suitable

Isn't that basically the Newton-Raphson algorithm without the iteration?

[FMDenaro]

Originally Posted by hunt_mat

Isn't that basically the Newton-Raphson algorithm without the iteration?

It depends on how you perform the iteration. The system becomes linear, you could suppose to develop a sort of Jacobi or Gauss-Seidel like iteration.
Of course, if you want a greater accuracy in time, the backward time-integration must be improved.


Your integration drives the solution to the rest condition?

[hunt_mat]

Originally Posted by FMDenaro

It depends on how you perform the iteration. The system becomes linear, you could suppose to develop a sort of Jacobi or Gauss-Seidel like iteration.
Of course, if you want a greater accuracy in time, the backward time-integration must be improved.


Your integration drives the solution to the rest condition?

Suppose my system of algebraic equations is defined as

F_j(u_{n+1},u_n)=0

Then I would define my linearization as:

F_j(u_n,u_n)+dF_j/du×du=0

This a now a linear system that I can solve using matlabs linear equation solver. I compute u_{n+1} by adding du to u_n.

[FMDenaro]

Originally Posted by hunt_mat

Suppose my system of algebraic equations is defined as

F_j(u_{n+1},u_n)=0

Then I would define my linearization as:

F_j(u_n,u_n)+dF_j/du×du=0

This a now a linear system that I can solve using matlabs linear equation solver. I compute u_{n+1} by adding du to u_n.

In your specific case:


u(i,n+1)=u(i,n)+dt*du/dt|(i,n)+...


du/dt(i,n)= -u(i,n)*du/dx|(i,n)+D* d2u/dx^2|(i,n)


but due to the poor time accuracy, I would start by simply setting the convective term in the form


u(i,n)*(u(i+1,n+1)-u(i,n+1))/dx




without further term.


But are you sure the velocity is always negative and the node i=1 remains at u=0?

[hunt_mat]

Originally Posted by FMDenaro

In your specific case:


u(i,n+1)=u(i,n)+dt*du/dt|(i,n)+...


du/dt(i,n)= -u(i,n)*du/dx|(i,n)+D* d2u/dx^2|(i,n)


but due to the poor time accuracy, I would start by simply setting the convective term in the form


u(i,n)*(u(i+1,n+1)-u(i,n+1))/dx




without further term.


But are you sure the velocity is always negative and the node i=1 remains at u=0?

The answers to your questions are yes to both. I am modelling a rod, standing on its head which is being heated from above, the density of the rod(that I have to include). The rod starts to heat up and starts to move from the free end due to the warming and as it heats up it moves to the end that is on the oven floor (that obviously doesn't move). As mass is conserved, the length of the rod has to decrease, so the velocity is towards the bottom, and therefore always negative.

[FMDenaro]

Originally Posted by hunt_mat

The answers to your questions are yes to both. I am modelling a rod, standing on its head which is being heated from above, the density of the rod(that I have to include). The rod starts to heat up and starts to move from the free end due to the warming and as it heats up it moves to the end that is on the oven floor (that obviously doesn't move). As mass is conserved, the length of the rod has to decrease, so the velocity is towards the bottom, and therefore always negative.

I am not sure to understand, I never worked on this kind of problem, you are saying that in your problem the u variable is the vertical velocity?


Check the textbook of LeVeque, I remember (if I am not wrong) that a similar problem is detailed.

[hunt_mat]

Originally Posted by FMDenaro

I am not sure to understand, I never worked on this kind of problem, you are saying that in your problem the u variable is the vertical velocity?


Check the textbook of LeVeque, I remember (if I am not wrong) that a similar problem is detailed.

I'm aware of the example but it's different to the problem I'm studying as it isn't the same, there are some very important differences such as a variable domain for one. Diffusion of the material is another.