# FEM for nonlinear first-order ODE

Currently I am trying to solve nonlinear Ricatti equation using FEM (Matlab language): $$\frac{d r(z)}{dz} = i k(z) r(z) + \frac{i k(z)}{2}(\epsilon(z) - 1)(1 + r(z))^2$$ $$z = [-h/2, h/2]$$ and $$r(-h/2)=0$$. $$k(z) = 2\cdot \pi\cdot n(z)/\lambda$$ $$n(z) = 1 + sin(\pi\cdot(\frac{z}{h} + \frac{1}{2}))$$ $$\varepsilon(z) = n^2(z)$$

First of all I write weak formulation: $$\int\frac{d r(z)}{dz} \cdot N(z) dz= \int i k(z) r(z)\cdot N(z)dz +\int \frac{i k(z)}{2}(\epsilon(z) - 1)(1 + r(z))^2\cdot N(z)dz$$

It should give me system of equations for undetermided coefficients (it is called weighted residuals method as I remember). At this point I have written code for FEM for this problem. It doesnt work correctly, as I am comparing results with ode45. I provide my code below:

clc
clear
h = 10;
n0 = 1;
n_elements = 5;
n_nodes = 2*n_elements+1;
L = h/n_elements;
z = -h/2:L/2:h/2;
tol = 10^(-5);
F = zeros(n_nodes,1);
F(1,1) = 0;
cnt = 0;
n = @(z) 1 + n0*sin(pi*(z/h +0.5));
lambda_range = 0.1*h:0.1*h:10*h;
lda = 10*h;
reflection = zeros(size(lambda_range));
tic
for j = 1:length(lambda_range)
lda = lambda_range(j);
while 1
F1 = assembly(F,n_elements,L,n,lda);
for i = 1:n_nodes
if abs(F(i)-F1(i))>tol
break;
end
end
F = F1;
cnt = cnt +1;
if cnt>1000
break
end
end
reflection(1,j) = (abs(F(end,1)))^2;
end
fprintf('Number of elements=%d\n',n_elements)
fprintf('Number of iterations=%d\n',cnt)

Plotting
plot(lambda_range,reflection,'b','Linewidth',3)
xlabel('\lambda')
ylabel('r(\lambda)')
title('Solution')

plot(z,abs(F).^2,'b','Linewidth',3)
xlabel('z')
ylabel('r(z)')
title('Solution')
toc

%%ODE45 check
tspan = [-0.5 0.5];
y0 = 0;
for j = 1:length(lambda_range)
lda = lambda_range(j);
[z,y] = ode45(@(z,y) 1i*2*pi*(1 + sin(pi*(z +0.5)))/lda*y+1i*2*pi*(1 + sin(pi*(z +0.5)))/lda*(((1 + sin(pi*(z +0.5))))^2-1)/2*(1+y)^2, tspan, y0);
reflection(1,j) = (abs(y(end,1)))^2;
end
plot(z,abs(y).^2,'b','Linewidth',3)
xlabel('z')
ylabel('r(z)')
title('Solution')
toc

plot(lambda_range,reflection,'b','Linewidth',3)
xlabel('z')
ylabel('r(z)')
title('Solution')
toc

%%Functions
function [F1] = assembly(F,n_elements,L,n,lda)
n_nodes = 2*n_elements+1;
K_table = [-1/2 -2/3 1/6; 2/3 0 -2/3; -1/6 2/3 1/2];
K = zeros(n_nodes,n_nodes);
R = zeros(n_nodes,1);
lmm = [];
for i =1:n_elements
j = (i-1)*2+1;
lmm = [lmm;[j,j+1,j+2]];
end

for i =1:n_elements
lm = lmm(i,:);

k11 = K_table-1i*(int1+int2) - 1i*(int_nl_1*F(lm(1)) + int_nl_2*F(lm(2)) + int_nl_3*F(lm(3)));
f1 = 1i*f11;

K(lm,lm) = K(lm,lm)+k11;
R(lm)=R(lm)+f1;
end

%%Boundary condition
K(1,:) = 0;
K(1,1) = 1;
R(1,1) = F(1);

%%Solution of a system
d = K\R;
F1 = d;
end

function [int1,int2, int_nl_1,int_nl_2,int_nl_3,f11, dop, k_table] = quadraticelement(L,i,n,lda)
gL = [-sqrt(3/5),0,sqrt(3/5)];
gW = [5/9, 8/9, 5/9];
k_table = zeros(3); int1 = zeros(3); int2 = zeros(3);int_nl_1 = zeros(3);int_nl_2 = zeros(3);int_nl_3 = zeros(3);f11 = sparse(3,1);

for k = 1:length(gW)
s = gL(k); w = gW(k);
N = [(s-1)*s/2, 1-s^2, (s+1)*s/2];
dns = [(2*s-1)/2, -2*s, (2*s+1)/2];
coord = [(i-1)*L; (i-1/2)*L; i*L];
z = L*s/2+(i-1/2)*L;
J = L/2;
dz = dns*(1/J);

f11 = f11 + J*w*N'*(2*pi*n(z)/lda*((n(z))^2 - 1)/2);
k_table = k_table + J*w*N.'*dz;
int1 = int1 + J*w*(N')*N*2*pi*n(z)/lda;
int2 = int2 + J*w*(N')*N*2*pi*n(z)/lda*((n(z))^2 - 1);
dop = [0, 0, 1];

int_nl_1 = int_nl_1 + J*w*(N')*N*N(1)*(2*pi*n(z)/lda*((n(z))^2 - 1)/2);
int_nl_2 = int_nl_2 + J*w*(N')*N*N(2)*(2*pi*n(z)/lda*((n(z))^2 - 1)/2);
int_nl_3 = int_nl_3 + J*w*(N')*N*N(3)*(2*pi*n(z)/lda*((n(z))^2 - 1)/2);
end
end

• For what values of $z$ are you trying to solve this equation? Nov 5, 2023 at 20:57
• @WolfgangBangerth z = [-h/2, h/2]. h can be 1 for example Nov 5, 2023 at 21:01
• Is this a boundary value problem or an initial value problem? If it is an IVP, why not use a standard ODE method (e.g., Runge-Kutta, BDF, etc.)? Nov 6, 2023 at 9:34
• It doesn't make sense to use the FEM for an initial-value problem. Nov 7, 2023 at 0:00
• Beside the fact that FEM is definitely not the right tool for this job: After discretization, you have $S \vec{u} = F(\vec{u})$, where $\vec{u}$ are the coefficients of the discrete solution $u_h$, matrix $S$ is the differential operator on the LHS and $F$ is the nonlinear RHS. You can use a Newton method in the discrete space to solve for $\vec{u}$, or try to solve by iteration. Nov 7, 2023 at 0:56