memory shortage: Solve Stiff ODEs

Dear matlab users.

I have memory shortage problem .

I use Solve Stiff ODEs https://jp.mathworks.com/help/matlab/math/solve-stiff-odes.html?lang=en

function brussode(N)
if nargin<1
   N = 100;
end
tspan = [0; 10];
y0 = [repmat(0.9,1,N); repmat(100200,1,N)];
options = odeset('Vectorized','on','JPattern',jpattern(N));
 
[t,y] = ode15s(@f,tspan,y0,options);
 
u = y(:,1:2:end);
x = (1:N)/(N+1);
p=y(:,2:2:end);
figure;
plot(x,u(end,:))
figure;
plot(x,p(end,:))

   function dydt = f(~,y)
a=0.0002/N;
b=3.55*10^-4;
c=2.39*10^-5;
k=3*10^-12; 

      dydt = zeros(2*N,size(y,2));      % preallocate dy/dt
      
      i = 1;
dydt(i,:) = 1/2/a^2/b*((k*(y(i,:).^3+y(i+2,:).^3)).*((y(i+3,:)-y(i+1,:))+ +19206*(1.417*(y(i+2,:)-y(i,:))-2.12*((1-y(i,:)).^2-(1-y(i+2,:)).^2)+1.263*((1-y(i,:)).^3-(1-y(i+2,:)).^3))).*y(i,:)-1.29*10^-9*a*b*2);
dydt(i+1,:) =1/2/a^2/c*((k*((1-y(i,:)).^3+(1-y(i+2,:)).^3).*(y(i+3,:)-y(i+1,:))).*y(i+3,:)+0.1./y(i+1,:) *a*c*2);   
      % Evaluate the 2 components of the function at all interior grid points.
  
i = 3:2:2*N-3;
dydt(i,:) = 1/2/a^2/b*((k*(y(i,:).^3+y(i+2,:).^3)).*((y(i+3,:)-y(i+1,:))+19206*(1.417*(y(i+2,:)-y(i,:))-2.12*((1-y(i,:)).^2-(1-y(i+2,:)).^2)+1.263*((1-y(i,:)).^3-(1-y(i+2,:)).^3))).*y(i,:)-(k*(y(i-2,:).^3+y(i,:).^3)).*((y(i+1,:)-y(i-1,:))+19206*(1.417*(y(i,:)-y(i-2,:))-2.12*((1-y(i-2,:)).^2-(1-y(i,:)).^2)+1.263*((1-y(i-2,:)).^3-(1-y(i,:)).^3))).*y(i-2,:));
dydt(i+1,:) =1/2/a^2/c*((k*((1-y(i,:)).^3+(1-y(i+2,:)).^3).*(y(i+3,:)-y(i+1,:))).*y(i+3,:)-(k*((1-y(i-2,:)).^3+(1-y(i,:)).^3)).*(y(i+1,:)-y(i-1,:)).*y(i+1,:));
 
      i = 2*N-1;
dydt(i,:) = 1/2/a^2/b*(1.29*10^-9*a*b*2-(k*(y(i-2,:).^3+y(i,:).^3)).*((y(i+1,:)-y(i-1,:))+19206*(1.417*(y(i,:)-y(i-2,:))-2.12*((1-y(i-2,:)).^2-(1-y(i,:)).^2)+1.263*((1-y(i-2,:)).^3-(1-y(i,:)).^3))).*y(i-2,:));
dydt(i+1,:) =1/2/a^2/c*(0.1./y(i+1,:) *a*c*2-(k*((1-y(i-2,:)).^3+(1-y(i,:)).^3)).*(y(i+1,:)-y(i-1,:)).*y(i+1,:));
    ;
   end
end
 
function S = jpattern(N)
% Jacobian sparsity pattern
B = ones(2*N,5);
B(2:2:2*N,2) = zeros(N,1);
B(1:2:2*N-1,4) = zeros(N,1);
S = spdiags(B,-2:2,2*N,2*N);
end

Please let me know what I need to change to calculate under 300GB.

Thank you for your help.

Answer 1

[Wan Ji]

你好！这个直接把

options = odeset('Vectorized','on','JPattern',jpattern(N));

改写成

options = odeset('Vectorized','on');

就可以了，因为你给出的ode方程的jpattern函数不能照搬别人的，所以干脆不用，这样计算也能快速得到想要的结果。

[颯太 小濱]

非常感谢你。 我立即试了一下，果然有效。

Thank you very much. I tried it immediately and it worked.