Facing problems in nonlinear system

%%%% Problem_01 %%%%
%u_t=u_xx+6u(1-u)%
%u(0,t)=(1+e^(-5t))^(-2)%
%u(1,t)=(1+e^(1-5t))^(-2)%
%u(x,0)=(1+e^x)^(-2)%

clc;
clear all;
format short

L = 1;       % Length of the rod
T = 0.05;       % Total time
Nx = 7;       % Number of spatial steps
Nt = 10;     % Number of time steps
alpha = 1;   % Thermal diffusivity
dx = L / Nx;   % Spatial step size
dt = T / Nt;   % Time step size
r = alpha * dt / dx^2
u = sym('u', [Nx+1,Nt+1]);

% Define the spatial grid
x = linspace(0, L, Nx+1);
x

% Set initial condition u(x,0)
u(:, 1) = vpa(0.0001679.*(x.^2-x)+0.00002215.*(1.5.*x.^3-1.5.*x),9);

% Set Dirichlet boundary conditions
u(1, :) = boundary_condition_x0(linspace(0, T, Nt+1)); % u(0,t)
u(end, :) = boundary_condition_xL(linspace(0, T, Nt+1)); % u(1,t)
u;

% Initialize the source term matrix
f = [];
R = [];
for n = 1:Nt
    for j = 1:Nx+1
        f_expr = (-10*(exp(x(j)-5*(n-0.5)*dt))/(1+exp(x(j)-5*(n-0.5)*dt))^3) -(2*exp(x(j)-5*(n-0.5)*dt)/(1+exp(x(j)-5*(n-0.5)*dt)^3)*((1-2*exp(x(j)-5*(n-0.5)*dt)))/(1+ ...
            exp(x(j)-5*(n-0.5)*dt)))-0.0003358-0.00019935*x(j) +6*((1/(1+exp(x(j)-5*(n-0.5)*dt))^2)-0.0001679*(x(j).^2-x(j))-0.00002215*((3/2)*(x(j).^3-x(j))))*(1-(1/(1+ ...
            exp(x(j)-5*(n-0.5)*dt))^2) +0.0001679*(x(j).^2-x(j))+0.00002215*((3/2)*(x(j).^3-x(j)))) + 6*0.5*(u(j,n)+u(j,n+1))*(1-0.5*(u(j,n)+u(j,n+1)));
        f{j,n} = f_expr;
    end
    for j = 2:Nx
        eq = (1-6*r)*u(j-1, n+1) + (10 + 12*r)*u(j, n+1) + (1 - 6*r)*u(j+1, n+1) == (1 + 6*r)*u(j-1, n)...
            + (10-12*r)*u(j, n) + (1 +6*r)*u(j+1, n) + dt*(f{j-1,n} +10*f{j,n} + f{j+1,n});
        eqs(n,j-1) = eq;
    end
    % disp("Equations before solving:");
    % disp(vpa(eqs(n, :), 6));
    vsol = vpasolve(eqs(n,:));
    R = struct2cell(vsol);
    for j = 2:Nx
        u(j,n+1) = min(abs(R{j-1}));
    end
end
vpa(u,9);

esol = @(x,t) (1+exp(x-5*t))^(-2);
exact_sol = [];
Compact_sol = [];
Error_Compact = [];
for n = 2:Nt+1
    for j = 2:Nx+1
        exact_sol(j,n) = esol((j-1)*dx,(n-1)*dt);
        Compact_sol(j,n) = u(j,n);
        Error_Compact(j,n) = abs(exact_sol(j,n)-Compact_sol(j,n));
    end
end
n = 11; % Choose any specific value of n (1 to 10)
j_values = 1:Nx+1;
u_values = u(j_values, n);
exact_val = exact_sol(j_values, n);
Compact_val = Compact_sol(j_values,n);
Compact_error_val = Error_Compact(j_values, n);

Table = table(u_values, exact_val,Compact_val,Compact_error_val, ...
    'VariableNames', {'E(i,j)', 'Exact_Solution','Compact_Solution','Compact_Error'})

function bc_x0 = boundary_condition_x0(t)% E(0,t)
    bc_x0 = zeros(size(t));
end

function bc_xL = boundary_condition_xL(t)% E(1,t)
    bc_xL = zeros(size(t));
end

[Torsten]

I don't understand your discretization in space and time. And it seems you implemented a different problem than the one you stated at the beginning of your code.

[Kashfi]

I use the compact finite difference method for both space and time discretizations, and I assure you that this discretization is correct. My goal is to improve the solution using this method. Initially, I solved the problem using the Weighted Residual Method, which is why the two approaches differ, but both are valid.

The main issue I am facing occurs when the number of spatial grid points, NxN_x, is 5 or less—the code produces a solution in this case. However, when I increase NxN_x, the code continues running indefinitely without producing a result.

[Torsten]

If you can choose which solver to use for your problem, I'd immediately choose "pdepe":

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

If you have to use your method, see the revised code below.

[Kashfi]

Actually I have to use my mentioned method, thats the challange for me.

Answer 1

[Torsten]

"pdepe" gives a different solution than your method. Maybe the nonlinear systems have multiple solutions.

%%%% Problem_01 %%%%
%u_t=u_xx+6u(1-u)%
%u(0,t)=(1+e^(-5t))^(-2)%
%u(1,t)=(1+e^(1-5t))^(-2)%
%u(x,0)=(1+e^x)^(-2)%

L = 1;         % Length of the rod
T = 0.05;      % Total time
Nx = 500;      % Number of spatial steps
Nt = 10;       % Number of time steps
alpha = 1;     % Thermal diffusivity
dx = L / Nx;   % Spatial step size
dt = T / Nt;   % Time step size
r = alpha * dt / dx^2;
x = linspace(0, L, Nx+1);
t = linspace(0, T, Nt+1);
% Set initial condition u(x,0)
u = zeros(Nx+1,Nt+1);
u(:,1) = 0.0001679*(x.^2-x)+0.00002215*(1.5*x.^3-1.5*x);
for n = 2:Nt
  u(:,n) = fsolve(@(U)fun(U,n-1,Nx,x,dx,t(n),dt,r,u(:,n-1)),u(:,n-1),optimset('Display','none'));  
end
figure(1)
plot(x,u)

uic = @(x)0.0001679*(x.^2-x)+0.00002215*(1.5*x.^3-1.5*x);
upde = @(x,t,u,dudx)deal(1,dudx,6*u*(1-u));
ubc = @(xl,ul,xr,ur,t)deal(ul,0,ur,0);
x = linspace(0,L,501);
t = linspace(0,T,11);
m = 0;
sol = pdepe(m,upde,uic,ubc,x,t);
u = sol(:,:,1);
figure(2)
plot(x,u)

●

function bc_x0 = boundary_condition_x0(t)% E(0,t)
    bc_x0 = zeros(size(t));
end

function bc_xL = boundary_condition_xL(t)% E(1,t)
    bc_xL = zeros(size(t));
end

function res = fun(U,n,Nx,x,dx,t,dt,r,Uold)
  f = zeros(Nx+1,1);
  res = zeros(Nx+1,1);
  for j = 1:Nx+1
    f(j) = (-10*(exp(x(j)-5*(n-0.5)*dt))/(1+exp(x(j)-5*(n-0.5)*dt))^3) -(2*exp(x(j)-5*(n-0.5)*dt)/(1+exp(x(j)-5*(n-0.5)*dt)^3)*((1-2*exp(x(j)-5*(n-0.5)*dt)))/(1+ ...
            exp(x(j)-5*(n-0.5)*dt)))-0.0003358-0.00019935*x(j) +6*((1/(1+exp(x(j)-5*(n-0.5)*dt))^2)-0.0001679*(x(j)^2-x(j))-0.00002215*((3/2)*(x(j)^3-x(j))))*(1-(1/(1+ ...
            exp(x(j)-5*(n-0.5)*dt))^2) +0.0001679*(x(j)^2-x(j))+0.00002215*((3/2)*(x(j)^3-x(j)))) + 6*0.5*(Uold(j)+U(j))*(1-0.5*(Uold(j)+U(j)));
  end
  res(1) = U(1) - boundary_condition_x0(t); 
  for j = 2:Nx
    res(j) = (1-6*r)*U(j-1) + (10 + 12*r)*U(j) + (1 - 6*r)*U(j+1) - ((1 + 6*r)*Uold(j-1)...
             + (10-12*r)*Uold(j) + (1 + 6*r)*Uold(j+1) + dt*(f(j-1) +10*f(j) + f(j+1)));
  end
  res(Nx+1) = U(Nx+1) - boundary_condition_xL(t);
end

[Kashfi]

Thank you for your assistance. The code now provides solutions; however, I require numerical solutions in the form of tables, as I initially specified in my first mentioned code. Additionally, I would like to observe the values of uuu after obtaining the solutions.

When setting Nx=500, the output matrix has dimensions 11×501, which is as expected. However, when changing Nx=10, the matrix dimensions remain unchanged. Could you please indicate where I should make modifications to ensure that the matrix size correctly reflects the value of Nx?

esol = @(x,t) (1+exp(x-5*t))^(-2);
exact_sol = [];
Compact_sol = [];
Error_Compact = [];
for n = 2:Nt+1
    for j = 2:Nx+1
        exact_sol(j,n) = esol((j-1)*dx,(n-1)*dt);
        Compact_sol(j,n) = u(j,n);
        Error_Compact(j,n) = abs(exact_sol(j,n)-Compact_sol(j,n));
    end
end

Unrecognized function or variable 'Nt'.

n = 11; % Choose any specific value of n (1 to 10)
j_values = 1:Nx+1;
u_values = u(j_values, n);
exact_val = exact_sol(j_values, n);
Compact_val = Compact_sol(j_values,n);
Compact_error_val = Error_Compact(j_values, n);
Table = table(u_values, exact_val,Compact_val,Compact_error_val, ...
    'VariableNames', {'E(i,j)', 'Exact_Solution','Compact_Solution','Compact_Error'})

Needed table in this form... thank you so much for your help.

[Torsten]

When setting Nx=500, the output matrix has dimensions 11×501, which is as expected. However, when changing Nx=10, the matrix dimensions remain unchanged. Could you please indicate where I should make modifications to ensure that the matrix size correctly reflects the value of Nx?

I cannot reproduce this problem. I just added your code part and everthing worked fine.

I might be mistaken, but I think the loops in your code part should start at 1 instead of 2:

for n = 1:Nt+1
    for j = 1:Nx+1
        exact_sol(j,n) = esol((j-1)*dx,(n-1)*dt);
        Compact_sol(j,n) = u(j,n);
        Error_Compact(j,n) = abs(exact_sol(j,n)-Compact_sol(j,n));
    end
end

instead of

for n = 2:Nt+1
    for j = 2:Nx+1
        exact_sol(j,n) = esol((j-1)*dx,(n-1)*dt);
        Compact_sol(j,n) = u(j,n);
        Error_Compact(j,n) = abs(exact_sol(j,n)-Compact_sol(j,n));
    end
end

In the code below, I further changed the loop to compute u from

for n = 2:Nt
  u(:,n) = fsolve(@(U)fun(U,n-1,Nx,x,dx,t(n),dt,r,u(:,n-1)),u(:,n-1),optimset('Display','none'));  
end

to

for n = 2:Nt+1
  u(:,n) = fsolve(@(U)fun(U,n-1,Nx,x,dx,t(n),dt,r,u(:,n-1)),u(:,n-1),optimset('Display','none'));  
end

in order to include the solution at t = T.

Could you include the correct PDE you are trying to solve with your code ? The Problem_01 description does not seem to be the one you try to code at the moment because initial and boundary conditions are different.

%%%% Problem_01 %%%%
%u_t=u_xx+6u(1-u)%
%u(0,t)=(1+e^(-5t))^(-2)%
%u(1,t)=(1+e^(1-5t))^(-2)%
%u(x,0)=(1+e^x)^(-2)%

L = 1;         % Length of the rod
T = 0.05;      % Total time
Nx = 500;      % Number of spatial steps
Nt = 10;       % Number of time steps
alpha = 1;     % Thermal diffusivity
dx = L / Nx;   % Spatial step size
dt = T / Nt;   % Time step size
r = alpha * dt / dx^2;
x = linspace(0, L, Nx+1);
t = linspace(0, T, Nt+1);
% Set initial condition u(x,0)
u = zeros(Nx+1,Nt+1);
u(:,1) = 0.0001679*(x.^2-x)+0.00002215*(1.5*x.^3-1.5*x);
for n = 2:Nt+1
  u(:,n) = fsolve(@(U)fun(U,n-1,Nx,x,dx,t(n),dt,r,u(:,n-1)),u(:,n-1),optimset('Display','none'));  
end
plot(x,u)

esol = @(x,t) (1+exp(x-5*t)).^(-2);
exact_sol = [];
Compact_sol = [];
Error_Compact = [];
for n = 1:Nt+1
    for j = 1:Nx+1
        exact_sol(j,n) = esol((j-1)*dx,(n-1)*dt);
        Compact_sol(j,n) = u(j,n);
        Error_Compact(j,n) = abs(exact_sol(j,n)-Compact_sol(j,n));
    end
end

n = 11; % Choose any specific value of n (1 to 10)
j_values = 1:Nx+1;
u_values = u(j_values, n);
exact_val = exact_sol(j_values, n);
Compact_val = Compact_sol(j_values,n);
Compact_error_val = Error_Compact(j_values, n);
Table = table(u_values, exact_val,Compact_val,Compact_error_val, ...
    'VariableNames', {'E(i,j)', 'Exact_Solution','Compact_Solution','Compact_Error'})

Table = 501x4 table
      E(i,j)      Exact_Solution    Compact_Solution    Compact_Error
    __________    ______________    ________________    _____________

    8.3323e-20       0.31604           8.3323e-20          0.31604   
    0.00019268       0.31549           0.00019268           0.3153   
    0.00038559       0.31494           0.00038559          0.31455   
    0.00057851       0.31438           0.00057851          0.31381   
    0.00077123       0.31383           0.00077123          0.31306   
    0.00096358       0.31328           0.00096358          0.31232   
     0.0011554       0.31273            0.0011554          0.31157   
     0.0013466       0.31218            0.0013466          0.31083   
      0.001537       0.31163             0.001537          0.31009   
     0.0017267       0.31108            0.0017267          0.30935   
     0.0019154       0.31053            0.0019154          0.30861   
     0.0021032       0.30998            0.0021032          0.30787   
       0.00229       0.30943              0.00229          0.30714   
     0.0024758       0.30888            0.0024758           0.3064   
     0.0026605       0.30833            0.0026605          0.30567   
     0.0028442       0.30778            0.0028442          0.30494   

function bc_x0 = boundary_condition_x0(t)% E(0,t)
    bc_x0 = zeros(size(t));
end

function bc_xL = boundary_condition_xL(t)% E(1,t)
    bc_xL = zeros(size(t));
end

function res = fun(U,n,Nx,x,dx,t,dt,r,Uold)
  f = zeros(Nx+1,1);
  res = zeros(Nx+1,1);
  for j = 1:Nx+1
    f(j) = (-10*(exp(x(j)-5*(n-0.5)*dt))/(1+exp(x(j)-5*(n-0.5)*dt))^3) -(2*exp(x(j)-5*(n-0.5)*dt)/(1+exp(x(j)-5*(n-0.5)*dt)^3)*((1-2*exp(x(j)-5*(n-0.5)*dt)))/(1+ ...
            exp(x(j)-5*(n-0.5)*dt)))-0.0003358-0.00019935*x(j) +6*((1/(1+exp(x(j)-5*(n-0.5)*dt))^2)-0.0001679*(x(j)^2-x(j))-0.00002215*((3/2)*(x(j)^3-x(j))))*(1-(1/(1+ ...
            exp(x(j)-5*(n-0.5)*dt))^2) +0.0001679*(x(j)^2-x(j))+0.00002215*((3/2)*(x(j)^3-x(j)))) + 6*0.5*(Uold(j)+U(j))*(1-0.5*(Uold(j)+U(j)));
  end
  res(1) = U(1) - boundary_condition_x0(t); 
  for j = 2:Nx
    res(j) = (1-6*r)*U(j-1) + (10 + 12*r)*U(j) + (1 - 6*r)*U(j+1) - ((1 + 6*r)*Uold(j-1)...
             + (10-12*r)*Uold(j) + (1 + 6*r)*Uold(j+1) + dt*(f(j-1) +10*f(j) + f(j+1)));
  end
  res(Nx+1) = U(Nx+1) - boundary_condition_xL(t);
end

[Kashfi]

The numerical solution appears to be running correctly. However, in the results table, only the values at the following specific points are required:u(0.1,0.05),u(0.2,0.05),…,u(1,0.05).

This corresponds to 10 discrete points. The value of Nx should be chosen accordingly.

The partial differential equation to be solved is given by:

with initial and boundary conditions, E(0,t) = E(1,t) = 0

is this code run for any nonlinear problems ?

[Kashfi]

"I think the loops in your code part should start at 1 instead of 2"

When i start at 1, it has error... "Index in position 2 is invalid. Array indices must be positive integers or logical values."

[Torsten]

The numerical solution appears to be running correctly. However, in the results table, only the values at the following specific points are required:u(0.1,0.05),u(0.2,0.05),…,u(1,0.05).

Then you should choose Nx such that 0.1, 0.2,..,0.9,1 are part of the grid, i.e. Nx should be a multiple of 10.

is this code run for any nonlinear problems ?

What do you mean by "any nonlinear problems" ?

It solves the Nx+1 nonlinear equations

U(1) - boundary_condition_x0(t) = 0; 
(1-6*r)*U(j-1) + (10 + 12*r)*U(j) + (1 - 6*r)*U(j+1) - ((1 + 6*r)*Uold(j-1)...
    + (10-12*r)*Uold(j) + (1 + 6*r)*Uold(j+1) + dt*(f(j-1) +10*f(j) + f(j+1))) = 0
(2 <= j <= Nx)
U(Nx+1) - boundary_condition_xL(t) = 0

for U in each time step.

When i start at 1, it has error... "Index in position 2 is invalid. Array indices must be positive integers or logical values."

I don't get this error. As you can see, it runs without problems with MATLAB online.

By the way: "pdepe" is not able to solve the problem you posted:

L = 1;
T = 0.05;
uic = @(x)0.001679*(x^2-x)+0.0002215*(1.5*x^3-1.5*x);
upde = @(x,t,u,dudx)deal(1,dudx,...
         6*u*(1-u)-...
         (1+exp(x-5*t))^(-2)-...
         0.001679*(x^2-x)-...
         0.0002215*(1.5*x^3-1.5*x)*2*exp(x-5*t)*(1-2*exp(x-5*t))/...
         ((1+exp(x-5*t))^3*(1+exp(x-5*t)))-...
         0.003358-0.0019935*x+...
         6*((1-exp(x-5*t))^(-2)-...
         0.001679*(x^2-x)-0.0002215*(1.5*x^3-1.5*x))*...
         (1-(1+exp(x-5*t))^(-2)+...
         0.001679*(x^2-x)+0.0002215*(1.5*x^3-1.5*x)));
ubc = @(xl,ul,xr,ur,t)deal(ul,0,ur,0);
x = linspace(0,L,1001);
t = linspace(0,T,11);
m = 0;
sol = pdepe(m,upde,uic,ubc,x,t);

Warning: Failure at t=1.000000e-04. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (2.168404e-19) at time t.

Warning: Time integration has failed. Solution is available at requested time points up to t=0.000000e+00.

u = sol(:,:,1);
plot(x,u)

●

[Kashfi]

It seems that pdepe is unable to solve the problem, but I was able to obtain solutions using a numerical approach with the code. I believe this is acceptable. Now, I need to test similar nonlinear problems using the same code to determine whether it provides solutions for those cases as well.

By the way, thank you very much for your effort—I truly appreciate it.