# Open to feedback for my code

1 view (last 30 days)
Lavorizia Vaughn on 13 Nov 2021
Commented: Steven Lord on 14 Nov 2021
hello, everyone, the code you see below is for a specific plot. my assignment has the following question:
Run Case 1 iwith the five stepsizes h = 0.05/2^(k1), k = 1, 2, 3, 4, and h = 0.001. Compute the value of θ2(t = 100) for each stepsize. Consider the last result the exact solution, and plot the four errors as a function of h in a loglog-plot. Estimate the order of convergence from the slope.
My report should contain the MATLAB commands used, the five values of θ2(t = 100), the loglog-plot, and the estimated slope.
Below is my code, and i was wondering if i could get some help on making it better. thank you.
th1=1;
th2=1;
w1=0;
w2=0;
hs(1)=[0.05];
hs(2)=[0.05/2];
hs(3)=[0.05/4];
hs(4)=[0.05/8];
hs(5)=[1/1000];
th2s=[]; %i feel like this could be better
for h=hs
y=[th1,th2,w1,w2];
N=100/h;
for i=1:N
k1=h*fpend(y);
k2=h*fpend(y + k1/2);
k3=h*fpend(y + k2/2);
k4=h*fpend(y + k3);
y = y + (k1 + 2*k2 + 2*k3 + k4) / 6;
end
th2s=[th2s,y(2)]; %i feel like this could be better
errors = abs(th2s(1:4) - th2s(5))
hs = hs(1:4);
loglog(hs,errors)
grid on
xlabel('Stepsize h')
ylabel('Error')
xlim([10^-3 10^-1])
ylim([10^-10 10^-5])
polyfit(log(hs), log(errors), 1)
Jan on 13 Nov 2021
It was a waste of time to post my answer there, if you have the above solution already.
You do not mean "h = 0.05/2(k−1)" but "h = 0.05 / 2^(k−1)".

Jan on 13 Nov 2021
[] is Matlab's operator for a concatenation. [0.05] concatenates the numnbre 0.05 with nothing, therefore this is a waste of time only. Nicer:
hs = [0.05, 0.05/2, 0.05/4, 0.05/8, 1/1000];
The outer loop is not closed by an end.
You do not only feel that "th2s=[th2s,y(2)];" is not optimal, but the editor displays a corresponding hint already: Do not let arrays grow iteratively. With a pre-allocation:
th1 = 1;
th2 = 1;
w1 = 0;
w2 = 0;
hs = [0.05, 0.05/2, 0.05/4, 0.05/8, 1/1000];
th2s = nan(size(hs));
for k = 1:numel(hs)
h = hs(k);
y = [th1,th2,w1,w2];
N = 100 / h;
for i = 1:N
k1 = h * fpend(y);
k2 = h * fpend(y + k1/2);
k3 = h * fpend(y + k2/2);
k4 = h * fpend(y + k3);
y = y + (k1 + 2 * k2 + 2 * k3 + k4) / 6;
end
th2s(k) = y(2);
end
errors = abs(th2s(1:4) - th2s(5))
hs = hs(1:4);
loglog(hs, errors)
grid on
xlabel('Stepsize h')
ylabel('Error')
xlim([1e-3, 1e-1])
ylim([1e-10, 1e-5])
polyfit(log(hs), log(errors), 1)
The run time does not matter in your case, but keep in mind that 10^-3 is an expensive poweroperation while 1e-3 is a cheap constant. This matters, if you call the code thousands of times.
Prefer spaces around the operators and separators between elements of vectors. This is less ambiguous:
a = [0 - 1i]
b = [0 -1i]
c = [0, -1i]
d = 2.*.2 % ???
##### 2 CommentsShowHide 1 older comment
Steven Lord on 14 Nov 2021
d = 2.*.2 % ???
d = 0.4000
In this case this would be written more clearly with the addition of one character. That extra character doesn't change its meaning but does make it easier for humans to parse.
d = 2.*0.2 % !
d = 0.4000

Image Analyst on 13 Nov 2021
Edited: Image Analyst on 13 Nov 2021
Y can be brought outside of the loop.
Here's a vectorized way to compute hs. Also note that for k = 1, (k-1) is zero, so hs is zero.
k = 1 : 4;
hs = (0.05 / 2) * (k-1);
hs(end+1) = 0.001
hs = 1×5
0 0.0250 0.0500 0.0750 0.0010

R2021a

### Community Treasure Hunt

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

Start Hunting!