MATLAB Answers

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?

9 views (last 30 days)
Dan on 6 Nov 2019
Commented: Dan on 7 Nov 2019
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);
i = 0;
while X(2,i+2)-X(2,i+1) > 0
i = i +1;
p = i;
while X(3,p+2)-X(3,p+1) > 0
p = p + 1;
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;
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;
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;
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)];
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')
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')
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')
grid on
hold on
title('Hermiteinterpoation - Kräfthatten')
hold on
legend('RK4','Linjen till spetsen', 'Hermiteinterpolation')
hold on
xlabel('Epsilon'); ylabel('Eta')


Show 1 older comment
Dan on 7 Nov 2019
I'm sorry for being inaccurate with the question.
The main question is how we can evaluate the accuracy of the results we get in this code?
Rik on 7 Nov 2019
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.

Sign in to comment.

Answers (1)

Bjorn Gustavsson
Bjorn Gustavsson on 7 Nov 2019
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:
This will now be a run in debug-mode, so you'll have a command-line prompt looking like this:
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.

Sign in to answer this question.


Translated by