memory shortage: Solve Stiff ODEs
2 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
颯太 小濱
el 18 de Ag. de 2021
Respondida: Wan Ji
el 18 de Ag. de 2021
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.
0 comentarios
Respuesta aceptada
Wan Ji
el 18 de Ag. de 2021
你好!这个直接把
options = odeset('Vectorized','on','JPattern',jpattern(N));
改写成
options = odeset('Vectorized','on');
就可以了,因为你给出的ode方程的jpattern函数不能照搬别人的,所以干脆不用,这样计算也能快速得到想要的结果。
Más respuestas (0)
Ver también
Categorías
Más información sobre Ordinary Differential Equations en Help Center y File Exchange.
Productos
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!