0
$\begingroup$

By method of lines I converted the PDE

u_t=u_{xx}, with the initial and boundary conditions

u(0,t)=u(0.1,t)
u(0,x)=sin(2*pi*x)

to a system of ODEs, as follows

function ut=pde1(t,u)
%
% Problem parameters
 % global ncall
  xl=0.0;
  xr=1.0;
  d=10;
  a=10;
%
% PDE
  n=length(u);
  h=((xr-xl)/(n-1));
  for i=1:n
    if(i==1)     ut(i)=0.0;
    elseif(i==n) ut(i)=0;
    else         ut(i)=d*(u(i+1)-2.0*u(i)+u(i-1))/h^2;
    end
  end
  ut=ut';
%

Now I am solving the ODEs ut by Runge-kutta method as follows:

%Initial condition
  n=500;
  for i=1:n
    u0(i)=sin((2*pi)*(i-1)/(n-1));
  end
  t0=0.0;
  tf=0.1;
  tout=linspace(t0,tf,n);
    y=zeros(n,length(tout));
 y(:,1)=u0';
    h=1/n;

    for i = 1 : length(tout)
        t=tout(i);
         k1 = pde_1(t,y(:,i));
         k2 = pde_1(t+h/2, y(:,i)+h*k1/2);
         k3 = pde_1(t+h/2, y(:,i)+h*k2/2);
         k4 = pde_1(t+h, y(:,i)+h*k3);
         y(:,i+1) = y(:,i) + h*(k1 + 2*k2 + 2*k3 + k4)/6;

    end

    y  
    plot(tout,y(:,1))

but the outputs are NAN! any help?

$\endgroup$
  • $\begingroup$ Why don't you use one of the Matlab built-in version of Runge-Kutta algorithm like ode45 ? $\endgroup$ – Jean Marie Jun 8 '17 at 2:24
  • $\begingroup$ @JeanMarie I am avoiding the build functions at this level for some a reason, I am considering some different types of runge kutta $\endgroup$ – Raed Jun 8 '17 at 4:36
  • $\begingroup$ Please correct the boundary conditions, order of variables, 0.1 vs. 1.0, is this a periodic b.c. or are both boundaries zero? If periodic, then your implementation is wrong. $\endgroup$ – LutzL Jun 8 '17 at 10:20
  • $\begingroup$ There are several problems with the code. As @LutzL pointed out, the problem statement is not consistent with the implementation. However, the reason behind your NaNs is that the CFL condition is not satisfied. For the diffusion equation you need the the time step $h$ to scale as $dx^2$ (where $dx$ is the spatial step), but in the code its simply set to $h = dx/10$. $\endgroup$ – ekkilop Jun 20 '17 at 15:50
  • $\begingroup$ @ekkilop you right, the stability condition should be satisfied $\endgroup$ – Raed Jun 23 '17 at 18:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.