How to plot a 2 variables function with integration and sum

3 vues (au cours des 30 derniers jours)
Hellema Ibrahim
Hellema Ibrahim le 1 Avr 2019
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

darova
darova le 1 Avr 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 commentaires
darova
darova le 2 Avr 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
Hellema Ibrahim
Hellema Ibrahim le 2 Avr 2019
Many thanks. Good day!

Connectez-vous pour commenter.

Plus de réponses (0)

Catégories

En savoir plus sur Numerical Integration and Differentiation dans Help Center et File Exchange

Produits


Version

R2016b

Community Treasure Hunt

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

Start Hunting!

Translated by