Open to feedback for my code

2 views (last 30 days)
Lavorizia Vaughn
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)
  1 Comment
Jan
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)".

Sign in to comment.

Accepted Answer

Jan
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 Comments
Steven Lord
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

Sign in to comment.

More Answers (1)

Image Analyst
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

Products


Release

R2021a

Community Treasure Hunt

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

Start Hunting!