- 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 view (last 30 days)
Show older comments
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 Comments
Answers (1)
Stephan
on 27 Nov 2018
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!