memory shortage: Solve Stiff ODEs
1 vue (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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 commentaires
Réponse acceptée
Wan Ji
le 18 Août 2021
你好!这个直接把
options = odeset('Vectorized','on','JPattern',jpattern(N));
改写成
options = odeset('Vectorized','on');
就可以了,因为你给出的ode方程的jpattern函数不能照搬别人的,所以干脆不用,这样计算也能快速得到想要的结果。
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Ordinary Differential Equations dans Help Center et File Exchange
Produits
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!