I am trying to solve a pair of coupled linear PDEs using backward Euler and Newton-Raphson
6 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
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 Kommentare
Torsten
am 24 Mai 2023
Bearbeitet: Torsten
am 24 Mai 2023
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.
Antworten (1)
Garmit Pant
am 22 Sep. 2023
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’:
I hope this helps!
Best Regards
Garmit
0 Kommentare
Siehe auch
Kategorien
Mehr zu Boundary Conditions finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!