solving a coupled pde by reducing it into coupled ode

Question

SDO il 16 Nov 2019

0
Link

Link diretto a questa domanda

https://it.mathworks.com/matlabcentral/answers/491438-solving-a-coupled-pde-by-reducing-it-into-coupled-ode

Risposto: SDO il 16 Nov 2019

Hi all

I am solving a highly nonlinear PDE y_t=[R(y,t)]_x by applying [R_(i+1/2) - R_(i-1/2)] / dx , reducing it to an ODE, and solving it with ode23tb, which works nicely.

However, when I try to do the same for a coupled system, y1=[R(y1,y2,t)]_x , y2=[Q(y1,y2,t)]_x it fails.

The sample code for the single equation is below, solving for periodic BC: ( y_t=[y^2y_xxx]_x )

Thanks in advance.

function sample
global N L
N = 51;L = 2*pi*wn;x = linspace(0,L,N);
st=10;ft=100;t = 0:st:ft;fin=ft/st;
y0 = 1-0.1*cos(x);
% sparsity matrix for the Jacobian
e = ones(N,1);S = spdiags([e e e e e], -2:2,N,N);
options = odeset('RelTol',1e-4,'AbsTol',1e-20, 'JPattern',S,'BDF','on');
[t,h] = ode23tb(@r,t,y0,options);
plot(x,h(fin,:),x,h(1,:),'--','LineWidth',2)
axis([0 L 0 2])
function yt = r(t,y)
global N L
yt0=size(N);
dx = L/(N-1);a2=-1/(2*dx^4);
r1=y.^2;
y(N+1) = y(2);y(N+2) = y(3);r1(N+1) = r1(2);
for i = 3:N
yt0(i)=a2*(r1(i+1) + r1(i))*(y(i+2)-3*y(i+1)+3*y(i)-y(i-1))-...
    a2*(r1(i) + r1(i-1))*(y(i+1)-3*y(i)+3*y(i-1)-y(i-2));
end
hm=y(N-1);hmm=y(N-2);r1m=r1(N-1);
i=1;
yt0(i) = (a2*(r1(i+1) + r1(i))*(y(i+2)-3*y(i+1)+3*y(i)-hm)-...
    a2*(r1(i) + r1m)*(y(i+1)-3*y(i)+3*hm-hmm));
i=2;
yt0(i) = (a2*(r1(i+1) + r1(i))*(y(i+2)-3*y(i+1)+3*y(i)-y(i-1))-...
    a2*(r1(i) + r1(i-1))*(y(i+1)-3*y(i)+3*y(i-1)-hm));
yt = yt0';