memory shortage: Solve Stiff ODEs

Question

颯太小濱 2021-8-18

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1435987-memory-shortage-solve-stiff-odes

回答： Wan Ji 2021-8-18

采纳的回答： Wan Ji

在 MATLAB Online 中打开

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