- you want to pass psivec to the function, but it is calculated inside the function.
- A is size 4096x4096 and wvec is size 64x1. Due to this A\wvec wont work.
What's wrong with my ode45 function?
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
I am using the ode45 to solve streamfunction and I compelete my functions and main script but always have wrong in ode 45 function, can anyone tell me where is it wrong, A,B,C are just matrix provided ( they are all correct). Thank you very much!
clear all,clc
% initial conditions
n=64;
tspan= 0:0.5:4;
v= 0.001;
L=20;
x2=linspace(-L/2,L/2,n+1);
xspan=x2(1:n);
dx = xspan(2)-xspan(1);
x=x2(1:n);
y=x;
w0=exp(-x.^2-((y.^2)/20));
wvec=reshape(w0,n,1);
% matrix A,B,C ( you can ignore this)
N=n*n; % total size of matrix
e0=zeros(N,1); % vector of zeros
e1=ones(N,1); % vector of ones
e2=e1; % copy the one vector
e4=e0; % copy the zero vector
for j=1:n
e2(n*j)=0; % overwrite every n^th value with zero
e4(n*j)=1; % overwirte every n^th value with one
end
e3(2:N,1)=e2(1:N-1,1); e3(1,1)=e2(N,1); % shift to correct
e5(2:N,1)=e4(1:N-1,1); e5(1,1)=e4(N,1); % positions
a=spdiags([e1 e1 e5 e2 2*e1 e3 e4 e1 e1], ...
[-(N-n) -n -n+1 -1 0 1 n-1 n (N-n)],N,N);
A=full(a)./((dx)^2);
b=spdiags([e1 -e1 e0 e1 -e1],[-n*(n-1) -n 0 n n*(n-1)],N,N);
B=full(b)./(2*dx);
c=spdiags([e5 -e2 e0 e3 -e4 ],[ -n+1 -1 0 1 n-1],N,N);
C=full(c)./(2*dx);
[t,wtsol]=ode45(@(t,wvec)rhseq(t,wvec,psivec,A,B,C,v),tspan,w0);
function wtsol=rhseq(t,wvec,psivec,A,B,C,v)
psivec=A\wvec;
wtsol=(C*psivec).*(B*wvec)-(B*psivec).*(C*wvec)+(v*A*wvec);
end
0 Kommentare
Antworten (1)
Stephan
am 27 Nov. 2018
Siehe auch
Kategorien
Mehr zu Ordinary Differential Equations finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!