FEM for nonlinear first-order ODE

3 次查看（过去 30 天）

显示更早的评论

Andrei 2023-11-6

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2043677-fem-for-nonlinear-first-order-ode

编辑： Andrei 2023-11-6

在 MATLAB Online 中打开

Currently I am trying to solve nonlinear Ricatti equation using FEM:

First of all I write weak formulation:

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

F = F1;

cnt = cnt +1;

if cnt>1000

break

end

reflection(1,j) = (abs(F(end,1)))^2;

end

fprintf('Number of elements=%d\n',n_elements)

Number of elements=5

fprintf('Number of iterations=%d\n',cnt)

Number of iterations=1100

%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

Elapsed time is 1.733690 seconds.

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

Elapsed time is 2.092498 seconds.

plot(lambda_range,reflection,'b','Linewidth',3)

xlabel('z')

ylabel('r(z)')

title('Solution')

toc

Elapsed time is 2.272228 seconds.

%%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,:);

[int1,int2, int_nl_1,int_nl_2,int_nl_3,f11] = quadraticelement(L,i,n,lda);

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

Where can be mistake, since for me code is correctly structured and methodology is also correct (or am I wrong?)?

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

Andrei 2023-11-6

Oh, I see. I tried to find an example and found it here (not a code and here linear problem is being observed). I tried to apply for my problem and extend it to nonlinear case with simple iterations.

Andrei 2023-11-6

编辑：Andrei 2023-11-6

Who wants to force you,

I study numerical optics this term. Basic task was to obtain solution using Runge-Kutta method of 2nd and 4th order, there was no problems. For "the best" students task was to obtain solution, but with FEM. Since previously I used to work in Fenics package on python (solved scattering from infinite cylinder) I was sure that using FEM will be easy))))))))

请先登录，再进行评论。

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