solving system of non linear PDEs

Question

Erm 2025-4-24

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2176608-solving-system-of-non-linear-pdes

编辑： Torsten 2025-4-25

在 MATLAB Online 中打开

Hello

I provide a solution to Allen cahn equation using spectral methods and Euler's formula.

The Allen Cahn Equation is

with BC

,

and BC

Now my question is if I have two systems of nonlinear partial differential equations (PDEs)., like

with BC

,

. How can I use the spectral method and Euler's formula to solve these problems?

Thank you. I truly appreciate any help

This script outlines the solution of the Allen-Cahn equation using spectral methods and Euler's formula:

% Allen-Cahn eq. u_t = u_xx + u - u^3, u(-1)=-1, u(1)=1
% (compare p6.m and p32.m).The challenge idea here is that the BC satisfy
% the IC. so Zeroing out the first and last rows of D2, This ensures that
% when you comput D2*(v - x) ,the first and last rows of this expression are
% always zero. So the time derivative at the boundary points is zero, i.e., those boundary values do not 
% change during time stepping, (  v + dt*(eps*D2*(v-x) + v - v.^3) = 0 at u(1) and u(end))
% Differentiation matrix and initial data:
N = 20; [D,x] = cheb(N); D2 = D^2; % use full-size matrix
D2([1 N+1],:) = zeros(2,N+1); % You're assigning values to two specific rows of the matrix D2 — the first and last.
                              % D2([1 N+1],:) means: "access rows 1 and N+1, and all columns. zeros(2,N+1) You're setting those rows
                              % to a matrix of zeros with size 2 by N+1
eps = 0.01; dt = min([.01,50*N^(-4)/eps]);
t = 0; v = .53*x + .47*sin(-1.5*pi*x);
% Solve PDE by Euler formula and plot results:
tmax = 100; tplot = 2; nplots = round(tmax/tplot);
plotgap = round(tplot/dt); dt = tplot/plotgap;
xx = -1:.025:1; vv = polyval(polyfit(x,v,N),xx);
plotdata = [vv; zeros(nplots,length(xx))]; tdata = t;
for i = 1:nplots
for n = 1:plotgap
t = t+dt; v = v + dt*(eps*D2*(v-x) + v - v.^3); % D2*(v-x): the function u(x,t)−x vanishes at the boundaries
end
vv = polyval(polyfit(x,v,N),xx);
plotdata(i+1,:) = vv; tdata = [tdata; t];
end
clf, subplot('position',[.1 .4 .8 .5])
mesh(xx,tdata,plotdata), grid on, axis([-1 1 0 tmax -1 1]),
view(-60,55), colormap([0 0 0]), xlabel x, ylabel t, zlabel u

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

Torsten 2025-4-25

编辑：Torsten 2025-4-25

在 MATLAB Online 中打开

I get substantially different results if I use the MATLAB standard solver "pdepe" for your problem. Maybe the cheb.m I added to your code is incorrect. If not, first take care that your code is able to solve only one PDE correctly.

% Allen-Cahn eq. u_t = u_xx + u - u^3, u(-1)=-1, u(1)=1

% (compare p6.m and p32.m).The challenge idea here is that the BC satisfy

% the IC. so Zeroing out the first and last rows of D2, This ensures that

% when you comput D2*(v - x) ,the first and last rows of this expression are

% always zero. So the time derivative at the boundary points is zero, i.e., those boundary values do not

% change during time stepping, ( v + dt*(eps*D2*(v-x) + v - v.^3) = 0 at u(1) and u(end))

% Differentiation matrix and initial data:

N = 20; [D,x] = cheb(N); D2 = D^2; % use full-size matrix

D2([1 N+1],:) = zeros(2,N+1); % You're assigning values to two specific rows of the matrix D2 — the first and last.

% D2([1 N+1],:) means: "access rows 1 and N+1, and all columns. zeros(2,N+1) You're setting those rows

% to a matrix of zeros with size 2 by N+1

eps = 0.01; dt = min([.01,50*N^(-4)/eps]);

t = 0; v = .53*x + .47*sin(-1.5*pi*x);

% Solve PDE by Euler formula and plot results:

tmax = 100; tplot = 2; nplots = round(tmax/tplot);

plotgap = round(tplot/dt); dt = tplot/plotgap;

xx = -1:.025:1; vv = polyval(polyfit(x,v,N),xx);

plotdata = [vv; zeros(nplots,length(xx))]; tdata = t;

for i = 1:nplots

for n = 1:plotgap

t = t+dt; v = v + dt*(eps*D2*(v-x) + v - v.^3); % D2*(v-x): the function u(x,t)−x vanishes at the boundaries

end

vv = polyval(polyfit(x,v,N),xx);

plotdata(i+1,:) = vv; tdata = [tdata; t];

end

%clf, subplot('position',[.1 .4 .8 .5])

%mesh(xx,tdata,plotdata), grid on, axis([-1 1 0 tmax -1 1]),

%view(-60,55), colormap([0 0 0]), xlabel x, ylabel t, zlabel u

pdefun = @(x,t,u,dudx)deal(1,dudx,u-u^3);

pdeic = @(x).53*x + .47*sin(-1.5*pi*x);

pdebc = @(xl,ul,xr,ur,t) deal(ul+1,0,ur-1,0);

m = 0;

sol = pdepe(m,pdefun,pdeic,pdebc,xx,tdata);

u = sol(:,:,1);

figure(1)

plot(xx,[plotdata(1,:);plotdata(2,:);plotdata(3,:);plotdata(end,:)])

figure(2)

plot(xx,[u(1,:);u(2,:);u(3,:);u(end,:)])

●

function [D,x] = cheb(N)

% CHEB compute D = differentiation matrix, x = Chebyshev grid, from

% Trefethen

if N==0, D=0; x=1; return, end

x = cos(pi*(0:N)/N)';

c = [2; ones(N-1,1); 2].*(-1).^(0:N)';

X = repmat(x,1,N+1);

dX = X-X';

D = (c*(1./c)')./(dX+eye(N+1));

D = D-diag(sum(D'));

end

Torsten 2025-4-25

编辑：Torsten 2025-4-25

You can replace the spatial derivatives by finite-difference approximations and then use the Euler method as time integrator (also for systems of PDEs) .

In case of the above PDE:

uxx(x_i ) ~ (u(x_(i+1))-2*u(x_i)+ u(x_(i-1)))/dx^2

for 2 <= i <= N-1,

thus you have to solve for x_new(x_i):

u_new(x_1) = -1

u_new(x_N) = 1

(u_new(x_i)-u_old(x_i))/dt = (u_old(x_(i+1))-2*u_old(x_i)+ u_old(x_(i-1)))/dx^2 + u_old(x_i) - u_old(x_i)^3

for 2 <= i <= N-1.

But usually the system above is stiff, and a very small time step dt is required for the Explicit Euler method in time to keep the solution stable.

Look up "Method-of-lines" for more information about the solution method. The advantage in using this method is that you can use a stiff MATLAB integrator like ODE15S for the time integration and that you only need to care about the spatial discretization.

Erm 2025-4-25

thank you. I appreciate your help.

I will read about MOF.

请先登录，再进行评论。

请先登录，再回答此问题。

solving system of non linear PDEs

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

回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

solving system of non linear PDEs

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

回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

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