Filter löschen
Filter löschen

newton raphson for nonlinear heat equation scheme

15 Ansichten (letzte 30 Tage)
Maryam Albuhayr
Maryam Albuhayr am 9 Sep. 2024
Kommentiert: Maryam Albuhayr am 9 Sep. 2024
I am trying to apply newton raphson method to linearize a scheme arise from non linear heat equation of the form
and the scheme is :
where are non linear term that was splitted to convex-concave function.
I have written this code but I couldn't get what is my mistake and how to make it work.
I appreciate any one can help me with that
close all
clear all
clc
%%
%PR parameters
R=8.314;
T=352.49;
% Critical constants for normal butane (c_6)
Tc = 507.44; % Critical temperature in Kelvin
Pc = 30.31; % Critical pressure in bar
w = 0.299; % Acentric factor
% MW=0.05812;
m=0.37464+1.54226*w-0.26992*w^2;
a=0.45724*R^2*Tc^2/Pc*(1+m*(1-sqrt(T/Tc)))^2;
b=0.07780*R*Tc/Pc;
%Influenceparameter
m1=10^(-16)/(1.2326+1.3757*w);
m2=10^(-16)/(0.9051+1.5410*w);
c=a*b^(2/3)*(m1*(1-T/Tc)+m2);
%% Parameters
L = 0.000000008; % Length of the domain
T = 0.1; % Total time
Nx = 5; % Number of spatial grid points
Nt = 10; % Number of time steps
% c = 1; % Thermal diffusivity constant
dx = L / (Nx - 1); % Spatial step size
dt = T / Nt; % Time step size
% Initial and boundary conditions
u0 = @(x) 7000-(7000/L)*x; % Initial condition
uL = 7000; % Dirichlet boundary condition at x=0
uR = 0; % Dirichlet boundary condition at x=L
% Nonlinear source term and its convex/concave splitting
f = @(u) R.*T*(log(u./(1-u.*b))+u.*b./(1-u.*b))-u.*a./(1+2.*u.*b-u.^2*b.^2)-a./(2.*sqrt(2).*b).*log((1+(1-sqrt(2)).*u.*b)./(1+(1+sqrt(2)).*u.*b)); % Example: cubic nonlinearity
fc = @(u) R.*T.*(log(u./(1-u.*b))+u.*b./(1-u.*b)); % Convex part
fd = @(u) f(u) - fc(u); % Concave part
fplot(f)
% Derivative of the convex part for Newton's method
dfc = @(u) R.*T*(1/(u.*(b.*u-1)).^2); % Derivative of the convex part
% Exact solution function (if available)
% For demonstration, assuming an exact solution function
% exact_solution = @(x, t) 1./(1+exp(x-t));
% Spatial grid
x = linspace(0, L, Nx);
% u=zeros(Nt,Nx);
% Initial condition
u = u0(x);
figure(1)
plot(x,u0(x))
grid on
title('Initial condition')
xlabel('x')
ylabel('u0')
%%
% Laplacian matrix (central difference)
e = ones(Nx,1);
A = spdiags([e -2*e e], -1:1, Nx, Nx) / dx^2;
% Apply boundary conditions to A
A(1,:) = 7000; A(1,1) = 7000; % Fix the boundary at x=0
A(end,:) = 0; A(end,end) = 0; % Fix the boundary at x=L 0;
%%
% Time integration loop
for n = 1:Nt
% Apply boundary conditions at each time step
u(1) = 0; % Boundary at x=0
u(Nx) = 0; % Boundary at x=L
% Initial guess for Newton's method (use the previous time step)
U = u;
% Newton iteration loop
for k = 1:100
% Calculate the residual F(U^k)
F_Uk = U - dt * c * A * U - dt * fc(U) - (u + dt * fd(u));
% Calculate the Jacobian J(U^k)
J_Uk = speye(Nx) - dt * c * A - dt * spdiags(dfc(U), 0, Nx, Nx);
% Solve for Newton update dU
dU = J_Uk \ (-F_Uk);
% Update U with the Newton step
U = U + dU;
% Check for convergence
if norm(dU, inf) < 1e-6
break;
end
end
% Update solution for the next time step
u = U;
% Optionally, plot the solution at selected times
if mod(n, Nt/10) == 0 % Plot every 10% of the total time steps
figure;
plot(x, u, '-o', 'DisplayName', 'Numerical Solution');
xlabel('x');
ylabel('u');
title(sprintf('Solution at time t = %.2f', n * dt));
grid on;
end
end
Error using *
Incorrect dimensions for matrix multiplication. Check that the number of columns in the first matrix matches the number of rows in the second matrix. To operate on each element of the matrix individually, use TIMES (.*) for elementwise multiplication.

Antworten (1)

Torsten
Torsten am 9 Sep. 2024
Bearbeitet: Torsten am 9 Sep. 2024
Why don't you use "pdepe" ? Ten lines of code maximum, and you are done.
To get rid of the error in the above code, replace
% Spatial grid
x = linspace(0, L, Nx);
by
% Spatial grid
x = linspace(0, L, Nx).';
But other problems follow.
Why do you set A as below ? Shouldn't this be u or U ?
A(1,:) = 7000; A(1,1) = 7000; % Fix the boundary at x=0
A(end,:) = 0; A(end,end) = 0; % Fix the boundary at x=L 0;
Are you sure about u(1) = 0 ? I thought you want to keep it at 7000.
% Apply boundary conditions at each time step
u(1) = 0; % Boundary at x=0
u(Nx) = 0; % Boundary at x=L

Produkte


Version

R2024a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by