Hello, we have received a Matlab code in school and the assignment is to explain what it means. I wonder if anybody could take us through the code and explain it?

Question

Dan 2019-11-6

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/489516-hello-we-have-received-a-matlab-code-in-school-and-the-assignment-is-to-explain-what-it-means-i-wo

评论： Dan 2019-11-7

clear; close all; clc;
 
R = 7.8;
a = 5.5;
b = 11; 
 
N = 8;
h = 2*pi/N;
u = 0:h:2*pi;
 
X = zeros(3,N); X(:,1) = [0;0;0];    
F = @(u,X) [(R - a*cos(u))/(sqrt(b^2 + (R - a*cos(u)).^2)); R*cos(X(1)); R*sin(X(1))];
 
for n = 1:N
    F1 = F(u(n), X(:,n));
    F2 = F(u(n) + h/2, X(:,n) + (h/2)*F1);
    F3 = F(u(n) + h/2, X(:,n) + (h/2)*F2);
    F4 = F(u(n) + h, X(:,n) + h*F3);
    X(:, n+1) = X(:,n) + (h/6)*(F1 + 2*F2 + 2*F3 + F4);
end
i = 0;
 
 
while X(2,i+2)-X(2,i+1) > 0
    i = i +1;
end
 
p = i;
while X(3,p+2)-X(3,p+1) > 0
    p = p + 1;
end
 
Cvek = zeros(4,N);
for i = 1:i
    HL = [X(3,i+1); X(3,i); X(1,i+1); X(1,i)];
    VL = [1 X(2,i+1) X(2,i+1).^2 X(2,i+1).^3;1 X(2,i) X(2,i).^2 X(2,i).^3;0 1 2*X(2,i+1) 3*X(2,i+1).^2;0 1 2*X(2,i) 3*X(2,i).^2];
    C = flip(VL\HL);
    Cvek(:,i) = C;
end
 
for i = 1:(p-i)
    HL = [X(3,(p+2)-i); X(3,(p+1)-i); -asin(sin(X(1,(p+2)-i))); -asin(sin(X(1,(p+1)-i)))];
    VL = [1 X(2,(p+2)-i) X(2,(p+2)-i).^2 X(2,(p+2)-i).^3;1 X(2,(p+1)-i) X(2,(p+1)-i).^2 X(2,(p+1)-i).^3;0 1 2*X(2,(p+2)-i) 3*X(2,(p+2)-i).^2;0 1 2*X(2,(p+1)-i) 3*X(2,(p+1)-i).^2];
    C = flip(VL\HL);
    Cvek(:,(p+1-i)) = C;  
end
 
for i = 1:(N-p)
    HL = [X(3,(N+2)-i); X(3,(N+1)-i); -asin(sin(X(1,(N+2)-i))); asin(sin(X(1,(N+1)-i)))];
    VL = [1 X(2,(N+2)-i) X(2,(N+2)-i).^2 X(2,(N+2)-i).^3;1 X(2,(N+1)-i) X(2,(N+1)-i).^2 X(2,(N+1)-i).^3;0 1 2*X(2,(N+2)-i) 3*X(2,(N+2)-i).^2;0 1 2*X(2,(N+1)-i) 3*X(2,(N+1)-i).^2];
    C = flip(VL\HL);
    Cvek(:,(N+1-i)) = C;  
end
 
P1 = [0; 0];
P2 = [0; sqrt(b^2 + (R - a)^2)];
xvek = [P1(1) P2(1) X(2,end)];
yvek = [P1(2) P2(2) X(3,end)];
 
plot(X(2,:),X(3,:),'*-')
hold on
plot(xvek, yvek,'*-')
hold on
for j = 1:i
    H = poly2sym(Cvek(:,j));
    fplot(H, [X(2,j) X(2,j+1)], 'blue')  
end
hold on
for j = 1:(p-i)
    H = flip(poly2sym(Cvek(:,(p+1-j))));
    fplot(H, [X(2,(p+2-j)) X(2,(p+1-j))], 'blue')
end
hold on
for j = 1:(N-p)
    H = flip(poly2sym(Cvek(:,(N+1-j))));
    fplot(H, [X(2,(N+2-j)) X(2,(N+1-j))], 'blue')
end
grid on
hold on
title('Hermiteinterpoation - Kräfthatten')
hold on
legend('RK4','Linjen till spetsen', 'Hermiteinterpolation')
hold on
xlabel('Epsilon'); ylabel('Eta')

4 个评论
显示 2更早的评论隐藏 2更早的评论

Rik 2019-11-7

Since the code does not have any comments, the variable names are not very descriptive, and you don't provide any context, your question is impossible to answer.

Have a read here and here. It will greatly improve your chances of getting an answer.

Dan 2019-11-7

okey, thank you!

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Bjorn Gustavsson 2019-11-7

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/489516-hello-we-have-received-a-matlab-code-in-school-and-the-assignment-is-to-explain-what-it-means-i-wo#answer_400342

At the matlab command-line prompt type:

dpstop in nameiofscript

where nameiofscript is the name of the script-file. Then you run the script at the matlab command-line:

nameiofscript

This will now be a run in debug-mode, so you'll have a command-line prompt looking like this:

K>>

and you can step throught the script line-by-line by pushing the button in the editor saying something like "step". Between steps you can inspect or plot different variables as you see fit. Don't modify them on the commandline since that might/would modify the workings of the script. This will allow you to describe what the script does.

HTH