How to plot a 2 variables function with integration and sum

3 views (last 30 days)
Hello, I would like to plot:
With :
I tried to break it down between the operations of integration and series sum.
% f1 = sin(n.*y.*pi);
% f2 = int(f1,y,0,L); or trapz(y, f1meshgrid,1);
% f3 = f2* exp(n.*t.*pi)*sin(n.*y.*pi);
% f4 = cumsum or symsum (f3,dimN)
L = 0.01;
NbPts = 10;
y = linspace(0,L,NbPts);
t = linspace(0,30,NbPts);
n = 0:1:NbPts;
f1 = @(n,y)sin(n.*y.*pi/L);
[N,Y] = meshgrid(n,y);
f1mtn = f1(N,Y);
f2 = trapz(n, f1mtn,2);
f3 = @(n,y,t) f2*exp(n.*t.*pi)*sin(n.*y.*pi);
[N2,Y2,T2] = meshgrid(n,y,t);
f3mtn = f3(N2,Y2,T2);
f4 = @(y,t) cumsum(f3mtn,3);
[Y3,T3] = meshgrid(y,t);
f4mty = f4(Y3,T3);
figure();
surfc(Y3,T,f4mty);
I obtain the error when I do f3mtn = f3(N2,Y2,T2) :
Error using * Inputs must be 2-D, or at least one input must be scalar. To compute elementwise TIMES, use TIMES (.*) instead.
But I have 3 variables and I don't see how I can reduce them. Any suggestion to code this? Many thanks.

Accepted Answer

darova
darova on 1 Apr 2019
clc, clear
L = 15;
NbPts = 50;
y = linspace(0,L,NbPts*3);
t = linspace(0,0.1,NbPts);
n = 1:2:NbPts;
[N1, Y1] = meshgrid(n,y);
cn = trapz(y, sin(pi*N1.*Y1/L) );
plot(cn)
[T2, Y2, Cn] = meshgrid(t,y,cn);
[~, ~, N2] = meshgrid(t,y,n);
k = N2*pi./L;
u = Cn.*exp(k.*T2).*sin(k.*Y2);
U = sum(u,3);
[T3,Y3] = meshgrid(t,y);
figure
surf(T3,Y3,U)
xlabel('Time axis')
ylabel('Y xis')
  3 Comments
darova
darova on 2 Apr 2019
We need to have all matrix the same dimension (3D). You can't multiply 1D x 3D or 2D x 3D
f1mtn = f1(N,Y); % matrix 2D
f2 = trapz(n, f1mtn,2); % Vector
f3 = @(n,y,t) f2*exp(n.*t.*pi)*sin(n.*y.*pi);
[N2,Y2,T2] = meshgrid(n,y,t);
% N2 - 3D matrix, Y2 - 3D matrix, T2 - 3D matrix
f3mtn = f3(N2,Y2,T2); % trying to multiply vector f2 by 3D
Maybe if you rewrite your code using loops it will be clearer. Two parts of the code below does the same
% [T2, Y2, Cn] = meshgrid(t,y,cn);
% [~, ~, N2] = meshgrid(t,y,n);
% k = N2*pi./L;
% u = Cn.*exp(k.*T2).*sin(k.*Y2);
u = zeros( length(t), length(y), length(cn) );
for i = 1:length(t)
for j = 1:length(y)
for k = 1:length(cn)
k0 = n(k)*pi/L;
u(i,j,k) = cn(k)*exp( k0*t(i) )*sin( k0*y(j) );
end
end
end
u = u';
Could you also clarify why we need 3 times more NbPts for y? - We dont. Size may be different

Sign in to comment.

More Answers (0)

Products


Release

R2016b

Community Treasure Hunt

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

Start Hunting!