Solving 2nd Order Couple PDE in Matlab

Working on recreation a mathematical model from a 2014 paper entitled "Viscoelasticity of Tau Proteins Leads to Strain Rate-Dependent Breaking of Microtubules during Axonal Stretch Injury: Predictions from a Mathematical Model." I’m trying to solve a coupled second order PDE using the pdepe function, but the equation I’m trying to solve has a coupled PDE in time. The pdepe anticipates the c value to be the multiplier for a single partial derivative, say U1 but I have partial U1 and Partial U2 w.r.t. time for both equations and I don’t think its accounting for them. Is there a way to incorporate the second partial derivative using this function, or how could I go about doing it. Here is my code.

x = [0 0.005 0.01 0.05 0.1 0.2 0.5 0.7 0.9 0.95 0.99 0.995 1]; t = [0 0.005 0.01 0.05 0.1 0.5 1 1.5 2]; m=0; sol = pdepe(m,@pde2014,@pde2014ic,@pde2014bc, x,t); u1 = sol (:,:,1); u2 = sol(:,:,2);

function [c,f,s] = pde2014(X,T,U,DUDX) eta = 0.35; %s %Lc = 0.22; %um Lc = 3.5714; L = 1; %um strate = 22.9; Z = diff(U2,t);

c = strate*[eta;eta]; f =[2*(Lc^2/L^2); 2*(Lc^2/L^2)] .*DUDX; s = [(U(2) - U(1))/L; (U(1) - U(2))/L]; end

function U0 = pde2014ic(X); U0 = [0;0];

function [pl,ql,pr,qr] = pde2014bc(Xl,Ul,Xr,Ur,T) pl = [Ul(1); 0]; ql =[0;1]; pr =[0;Ur(2)-T]; qr = [1;0];

[Torsten]

Please include the equations in a mathematical notation.

Best wishes

Torsten.

[Sabahdo]

The equations are:

-2*(L_c^2 / L^2)* (d^2U1(X,T)/dX^2) = (n/L)*(dU2(X,T)/dt - dU1(X,T)/dt) + (U2(X,T) - U1(X,T))/L

and

2*(L_c^2 / L^2)* (d^2U2(X,T)/dX^2) = (n/L)*(dU2(X,T)/dt - dU1(X,T)/dt) + (U2(X,T) - U1(X,T))/L

and I know all the constants, IC, and BCs.

[Torsten]

If you substitute V=U2-U1 and add the two equations from above, you get a single equation in V:

n/L*dV/dt + v/L = (L_c^2/L^2)*d^2V/dX^2

Best wishes

Torsten.

Answer 1

[Faycal bouzouidja]

As 𝑐 is a functionsuch that 𝑐 = 𝑐(𝑥, 𝑡); 𝐷, 𝑈 and 𝑘 are real constants.

𝜕𝑐

𝜕2𝑐

𝜕𝑐

𝜕𝑡 = 𝐷𝜕𝑥2 − 𝑈 𝜕𝑥 − 𝑘𝑐

0 ≤ 𝑥 ≤ 10

Initial condition  𝑐(𝑥, 0) = 0

∆𝑡=0.005s, ∆𝑥=0.1 m, 𝐷 = 1, 𝑈 = 5, 𝑘 = 0.5

Write a MATLAB code to solve the given boundary value problem numerically using an explicitfinite difference method, for the boundary conditions given below:

• Constant boundary conditions  𝑐(0, 𝑡) =100 and 𝑐(10, 𝑡) = 0

Plot solutions at t = 0.25, 0.5, 0.75 and 1 s.

• Derivative boundarycondition at the first grid point:

𝜕𝑐

= 1

𝜕𝑥x=O

𝑐(10, 𝑡) =0

Plot solutions at t = 2, 4, 6 and 8 s.

• Derivative boundarycondition at the last grid point:

𝑐(0, 𝑡) =100

𝜕𝑐

= 1

𝜕𝑥x=1O

Plot solutions at t = 0.5, 1, 1.5 and 2 s.

Use first-order central differencing for the derivative boundary