Open to feedback for my code

1 次查看(过去 30 天)
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 个评论
Jan
Jan 2021-11-13
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
Jan 2021-11-13
[] 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 个评论
Steven Lord
Steven Lord 2021-11-14
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

请先登录,再进行评论。

更多回答(1 个)

Image Analyst
Image Analyst 2021-11-13
编辑:Image Analyst 2021-11-13
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

类别

Help CenterFile Exchange 中查找有关 Logical 的更多信息

产品


版本

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by