Hello everyone I hope you are all doing well. Please I have a concern about the coupling of an ODE and an EDP.
9 views (last 30 days)
I would like to couple a second order differential equation (ODE) with the heat equation (PDE). For this, I used the RK4 method for the ODE and the Crank-Nicholson methode for the EDP. I don't know why my code does not give me a good results. Please give me a hand. I have attached two files to this request: one which presents my code and the other which shows the equation system
% Matlab code : The dynamic of slider using the heat equation
% I would like to use the Cranck Nicholson methode to solve the heat equation and
% incorporate it simultaneously into the second order of ODE
% Parameters to define the dynamic of slider.
Kr = 10; m = 10; F0 = 1; nu0 = 10; Vp = 10.^(-11); gama = 0.8; beta = 0.2; uc = 0.5;
% Parameters to define the heat equation.
rho = 2700; cs = 1000; Kc = 2.7;
% Initial temperature of heat
Ti = 150;
% boundary condition of heat equation.
Tc = 1050; Th = 1150;
% final lenght and time.
L = 2; % thickness of conduction chamber.
tf =80; % Final time.
% Parameters needed to solve the heat equation.
maxk = 100; % Number of time steps
dt = tf/maxk; % lenght on time steps
n = 100; % Number of space steps
dx = L/n; % lenght on space steps
% Definition of time
t = 0:dt:maxk.*dt;
% parameters used to make the reduction of another parameters.
H = L./2; % the half thickness of conduction chamber
Kd = Kc./(rho.*cs); % thermal diffusivity
alpha = (Kd.*dt)./((dx).^(2)); % reduction parameter
alpha1 = (alpha*dx*dx)./(Kc*H); % reduction parameter
% % Previsional matrice to incorporate data.
u = zeros(1,maxk+1);
v = zeros(1,maxk+1);
T = zeros(n+1,maxk+1);
% Defining the Matrices M_left :(ML)
aal(1:n-2) = -alpha;
bbl(1:n-1) = 2 + 2.*alpha;
ccl(1:n-2) = -alpha;
ML = diag(bbl,0)+diag(aal,-1)+diag(ccl,1);
% Defining the Matrices M_rigth :(ML)
aar(1:n-2) = alpha;
bbr(1:n-1) = 2 - 2.*alpha;
ccr(1:n-2) = alpha;
MR = diag(bbr,0)+diag(aar,-1)+diag(ccr,1);
% initial slip, velocity and enery.
u(1) = 0; % slip.
v(1) = 0; % velocity.
Tq(1) = Th; % temperature
% Definition of the three functions using by Rk4.
pascak1 = @(t,u,v) v;
pascak2 = @(t,u,v) (-1./m).*((Kr-gama).*u + nu0.*exp(beta*((Th-Tq)./(Th-Tc))).*v + F0 - Kr.*Vp.*t);
for i = 1:n+1
% Initial temperature of the fault.
T(i,1) = Ti;
% Boundary condition
for k=2:maxk + 1 % Time Loop
T(1,k) = Th;
T(n+1,k) = Tc;
% time(k) = (k-1)*dt;
L1(k) = dt.*pascak1(t(k-1),u(k-1),v(k-1));
K1(k) = dt.*pascak2(t(k-1),u(k-1),v(k-1));
L2(k) = dt.*pascak1(t(k-1) + dt/2, u(k-1) + L1(k)./2 ,v(k-1) + K1(k)./2);
K2(k) = dt.*pascak2(t(k-1) + dt/2, u(k-1) + L1(k)./2 ,v(k-1) + K1(k)./2);
L3(k) = dt.*pascak1(t(k-1) + dt/2, u(k-1) + L2(k)./2 ,v(k-1) + K2(k)./2);
K3(k) = dt.*pascak2(t(k-1) + dt/2, u(k-1) + L2(k)./2 ,v(k-1) + K2(k)./2);
L4(k) = dt.*pascak1(t(k-1) + dt, u(k-1) + L3(k) ,v(k-1) + K3(k));
K4(k) = dt.*pascak2(t(k-1) + dt, u(k-1) + L3(k) ,v(k-1) + K3(k));
% Implementation of Crank-Nicholson method
T(2:n,k)= inv(ML)*MR*TT + alpha1*F0*inv(ML)*v(k-1)*ones(99,1) - alpha1*gama*inv(ML)*u(k-1)*v(k-1)*ones(99,1);
% Because T is two dimenssions I extract Tq in one dimenssion at a specific space and
% incorporate it in the ODE equation.
Tq = T(50,1:maxk+1);
u(k) = u(k-1) + (1/6).*(L1(k) + 2*L2(k) + 2*L3(k) + L4(k));
v(k) = v(k-1) + (1/6).*(K1(k) + 2*K2(k) + 2*K3(k) + K4(k));
t(k) = t(k-1) + dt;
Bill Greene on 28 Sep 2022
I looked at your mathematical description of this problem.Since T is a function of z, V and u are also functions of z. So you really are solving three PDE. I see that in your code, you use T at the mid-point of the region in the second equation to convert it to an ODE. Are you sure this doesn't introduce significant errors?