What's wrong with my ode45 function?

1 view (last 30 days)
Yuxing Zhang
Yuxing Zhang on 27 Nov 2018
Commented: Yuxing Zhang on 27 Nov 2018
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

Answers (1)

Stephan
Stephan on 27 Nov 2018
  1. you want to pass psivec to the function, but it is calculated inside the function.
  2. A is size 4096x4096 and wvec is size 64x1. Due to this A\wvec wont work.

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!