Analytical solution for a difficult nonlinear PDE

Question

Is it possible to compute the analytical solution for this nonlinear pde? It doesn't seems to work with Sympy but it doesn't i can do it with it. The point is to prove that the convergeance order of the code below is indeed 2 (since it's using a newton method). Since i didn't manage to get the analytical solution at first, i tried to test a function that verify the boundary conditions (u(1) = 1 and u'(0) = 0) but it diverges. The function was uex = -cos(pi*X)

close all
clear all
clc
clear;
N=51; L=1;
dx=L/(N-1); X=[0:dx:L];
K0=0.01; sigma=1; beta=300;
I=2:N-1;

% uex = -cos(pi*X)
% Q = beta*4*(-cos(pi*X)-1)-0.01*(pi^2)*cos(pi*X)
delta=0.2; 
Q=(X<delta)*beta;

K=inline('t.^2'); 

Un=ones(1,N); dUn=zeros(1,N); Fn=zeros(1,N);

nitmax=1000; epsilon=1.0e-5;
for it=1:nitmax
    Kn12=0.5*(K(Un(1:N-1))+K(Un(2:N))); 
    Fn(I)=K0/dx^2*(Kn12(I-1).*(Un(I-1)-Un(I))+Kn12(I).*(Un(I+1)-Un(I)))-sigma*(Un(I).^4-1.0)+Q(I);
    Fn(1)=K0/dx^2*(Kn12(1)*(Un(2)-Un(1))+Kn12(1)*(Un(2)-Un(1)))-sigma*(Un(1).^4-1.0)+Q(1);
    Fn(N)=0;
    A=[K0/dx^2*(2*Kn12(1))+4*sigma*Un(1)^3,K0/dx^2*(Kn12(I-1)+Kn12(I))+4*sigma*Un(I).^3,1];
    B=[-K0/dx^2*(2*Kn12(1)), -K0/dx^2*(Kn12(I)),0];
    C=[0,-K0/dx^2*(Kn12(I-1)),0]; 
    J = full(gallery('tridiag', B, [A 0], C));
    J(end-1, 1) = J(end-1, end); 
    J(1, end-1) = J(end, end-1); 
    J = J(1:end-1, 1:end-1); 
    dUn=J\Fn.';
    Err=norm(Fn)/sqrt(N);
    Un=Un+dUn;  
    if (Err<epsilon) break; end;
end;
Un(:,1)
plot(Un(:,1))

Answer 1

Here is a scheme that might be useful, to reduce the problem to first order equations.

Take the case where $\kappa(u) = \kappa_0 u^2$ I will write primes for $x$ derivatives since this is not a PDE. Since $$ u^2 u' = \tfrac{1}{3}(u^3)' $$ your equation takes the form $$ c_1(u^3)''+\sigma u^4 = q(x) $$ where $q$ takes just two values on two subintervals of $[0,1]$, and for simplicity I won't track all the constants.

Substitute $y = u^3$, and multiply the equation by $y'$ to get $$ c_2 y'' y' +\sigma y^{4/3} y' = q(x)y'. $$ This is, on the two subintervals, $$ (c_3y'^2+c_4 y^{7/3})' = (q(x)y)'. $$ Integrate to get $$ c_3y'^2+c_4 y^{7/3} = q(x)y + c_5 $$ on the first subinterval and $$ c_3y'^2+c_4 y^{7/3} = q(x)y + c_6 $$ on the second.

These first order equations can possibly be solved numerically the first to match the boundary condition $u_x(0) = 0$ leaving $c_5$ as a parameter, the second for $u(1) = 1$, and finally chosing $c_5$ and $c_6$ to attempt to make $u$ continuous or maybe smooth at $x = 0.2$.