Facing problems in nonlinear system

Question

Kashfi 2025-3-28

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2175740-facing-problems-in-nonlinear-system

移动：Torsten 2025-3-30

%%%% 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

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

Torsten 2025-3-28

编辑：Torsten 2025-3-28

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 2025-3-29

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

请先登录，再进行评论。

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

Answer 1

Torsten 2025-3-28

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2175740-facing-problems-in-nonlinear-system#answer_1562787

编辑：Torsten 2025-3-28

在 MATLAB Online 中打开

"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

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

Kashfi 2025-3-29

移动：Torsten 2025-3-29

在 MATLAB Online 中打开

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 2025-3-29

编辑：Torsten 2025-3-29

在 MATLAB Online 中打开

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

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

Torsten 2025-3-29

编辑：Torsten 2025-3-29

在 MATLAB Online 中打开

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 2025-3-30

移动：Torsten 2025-3-30

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.

请先登录，再进行评论。

Facing problems in nonlinear system

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

采纳的回答

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

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

Facing problems in nonlinear system

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

采纳的回答

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

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

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

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