How to plot a 2 variables function with integration and sum

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.

 Réponse acceptée

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 commentaires

Thank you very much Darova. Would you mind to explain why the magical 2 lines:
[T2, Y2, Cn] = meshgrid(t,y,cn);
[~, ~, N2] = meshgrid(t,y,n);
I don't understand why we have to create a new vector Cn and new vector N2. Could you also clarify why we need 3 times more NbPts for y?
Thank you in advance.
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
Many thanks. Good day!

Connectez-vous pour commenter.

Plus de réponses (0)

Catégories

Produits

Version

R2016b

Community Treasure Hunt

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

Start Hunting!

Translated by