ODE45 Method of lines For two Coupled PDE

Hello all;

I am trying to solve numerically the following coupled differential equations using the method of lines:

du(t,y)/dt= Gr*Pr* Theta(t,y) + Pr * d2u(t,y)/dt^2

d(Theta(t,y))/dt= d2(Theta(t,y)/dy^2 + R*[(Theta(t,y)+Beta)^4 –Beta^4]

where Gr, Pr, and Beta are all constant:

First I discretized the equation using 2nd order central difference in space (y), then I solve for the transient using the ODE45 of Matlab. Discretization include applying the boundary conditions. Somehow the program stopped the integration over time after some time interval?? I don’t what is wrong. Obviously when someone uses the method of lines with two variables the interval of solution for the space has to be doubled i.e. the vector y in the script includes y(1:N+1) and Theta y(N+2:2N+3).

Can anyone of you guys help me of what is wrong?

Here are the codes

close all;

clear all;

n=10;

Fs =1; %sampling rate

Ts = 1/Fs; %sampling time interval

Per=n;dx=1/n;

alfa=2e-3;

%aa=alfa/dx^2;

aa=alfa;

tspan =0:Fs:Per; %sampling period

x=0:1/n:1;

%y0=zeros(1,2*(n+1));

y0(1,1:n+1)=0;

y0(1,n+2:2*n+3)=0;

[t,y] = ode45(@(t,y)rhs(t,aa,n,x,y),tspan,y0);

yy=y(:,1:n+1);yT=y(:,n+3:2*n+3);

size(y)

size(yy)

size(yT)

size(x)

for i=1:n+1

figure(1);

plot(x,yy(i,1:n+1),'k');hold on

%figure(2);

%plot(x,yT(i,1:n+1),'k');hold on

end

*****************************************

function ydot=rhs(t,aa,n,x,y)

dx=1/n;

u_bcl=1;u_bcr=0; %boundary condtion of u

T_bcl=1;T_bcr=0; %boundary condtion of Theta

beta=10;

R=1;

Pr=0.7;

Gr=10;

ydot=zeros(2*n+3,1);

%y(1)=0;y(n+1)=1;

y(n+2)=0;

%y(2*n+2)=1;

ydot(1) = (Pr/dx^2)*(y(2)-y(1)+u_bcl)+Gr*Pr*T_bcl; %first equation include BC for u

ydot(n+3)=(1/dx^2)*(y(n+4)-2*y(n+3)+T_bcl)+R*((y(n+3)+beta)^4-beta^4); %first equation include BC for Theta

for i =2:n

ydot(i) = (Pr/dx^2)*(y(i+1)-2*y(i)+y(i-1))+Gr*Pr*y(i+n+2);

end

for j=n+4:2*n+2

ydot(j) = (1/dx^2)*(y(j+1)-2*y(j)+y(j-1))+R*((y(j)+beta)^4-beta^4);

end

ydot(n+1) = (Pr/dx^2)*(u_bcr-2*y(n) + y(n-1))+Gr*Pr*T_bcr;

ydot(2*n+3)=(1/dx^2)*(T_bcr-2*y(2*n+3)+y(2*n+2))+R*((y(2*n+3)+beta)^4-beta^4);

[Torsten]

Use MATLAB's "pdepe".