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

Question

Matthew Hunt 2023-5-23

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1971544-i-am-trying-to-solve-a-pair-of-coupled-linear-pdes-using-backward-euler-and-newton-raphson

回答： Garmit Pant 2023-9-22

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.

4 个评论
显示 2更早的评论隐藏 2更早的评论

Torsten 2023-5-24

编辑：Torsten 2023-5-24

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 2023-5-25

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 2023-9-22

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1971544-i-am-trying-to-solve-a-pair-of-coupled-linear-pdes-using-backward-euler-and-newton-raphson#answer_1315717

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:

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.
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)
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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

4 个评论
显示 2更早的评论隐藏 2更早的评论

回答（1 个）

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

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

4 个评论 显示 2更早的评论隐藏 2更早的评论

回答（1 个）

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

4 个评论
显示 2更早的评论隐藏 2更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论