please give me your suggestions. thanks
Hello everyone I hope you are all doing well. Please I have a concern about the coupling of an ODE and an EDP.
2 views (last 30 days)
Show older comments
SARSKOLIN FOSSO on 26 Sep 2022
Commented: SARSKOLIN FOSSO on 29 Sep 2022
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?
Find more on Mathematics in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!Start Hunting!