How to plot a 2 variables function with integration and sum

6 次查看(过去 30 天)
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.

采纳的回答

darova
darova 2019-4-1
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 个评论
darova
darova 2019-4-2
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

请先登录,再进行评论。

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Numerical Integration and Differentiation 的更多信息

产品


版本

R2016b

Community Treasure Hunt

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

Start Hunting!

Translated by