Open to feedback for my code

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^(k−1), 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]

See your other question: https://www.mathworks.com/matlabcentral/answers/1585584-help-with-a-loop-that-computes-at-each-step

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)".

Answer 1

[Jan]

[] 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  % ???

[Lavorizia Vaughn]

ty

[Steven Lord]

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

Answer 2

[Image Analyst]

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