euler method first order

13 ビュー (過去 30 日間)
dulanga
dulanga 2019 年 4 月 1 日
回答済み: dulanga 2019 年 4 月 1 日
This is my code is this correct?
y0 = -1; % Initial Condition
h = 0.1;
t = 0:h:100;
N=length(t)
y_euler = zeros(size(t));
% Initial condition gives solution at t=0.
y_euler(1) = y0;
% Solving the equation via Euler's method
for i=1:(length(t)-1)
k1 = sqrt(t(i)*(y_euler(i)^2))-sqrt(t(i));
y_euler(i+1) = y_euler(i) + h*k1;
end
N50=(N-1)/2

採用された回答

John D'Errico
John D'Errico 2019 年 4 月 1 日
編集済み: John D'Errico 2019 年 4 月 1 日
Um, no? Not correct. Not terrible. Even pretty close. But not correct. You got a lot of it right. You preallocated y_euler correctly.
You should learn a few things, but just mainly minor points:
Get used to using semi-colons at the end of lines. You missed one when you created N, and another for N50.
You compute N, but then don't bother to use it, since you then use length(t) again.
White space is never a bad thing in your code. There is no charge for blank lines.
You create N50, but then never seem to use it. Be careful if N was even, in case N50 is intended to be later used as an index.
More comments are always better than just a few. They help you to read your code, like a good book. That will be important to you when you need to debug your code one day.
You might have chosen a more descriptive name for N. Good code is easy to read. It helps you to remember what each variable means in context. So a name like Nt might be more suggestive that Nt is the length of t. Similarly, I might have used a name for h as delta_t, or some such thing. But h is also a good choice, since that is commonly seen in the formulas for Euler's method as the step. The important thing is to use mnemonic names, so when you are reading through a long code in the future to debug, they remind you what each variable means.
Instead of hard coding t = 0:h:100, you might have written it more descriptively as:
t0 = 0;
tend = 100;
t = t0:h:tend;
Be careful though when you use colon, as unless (tend - t0)/h is an integer, the last step in the list created by colon may not indeed be tend. For example
0:0.1:pi
will not leave the last element as exactly pi, but as only 3.1. At some point in your coding future, that will become important.
Next, one slightly more significant point:
You hard-coded the function computation in your loop. This is a bad idea, because then it makes it hard to change your code to allow a different function at some point. Get used to using and creating functions!!!! Function handles are easy to create on the fly. Also, see that when you need to write code for a more complex method like Runge-Kutta soon, having dydt as a separate function handle will make that code greatly easier to read AND write.
Finally, you made one major mistake:
The function you use for k1, is significantly different from that which your assignment indicates. You wrote:
sqrt(t(i)*(y_euler(i)^2))-sqrt(t(i))
Is that really supposed to be
sqrt( t + y^2) - sqrt(t)
Well, I suppose it is close. You replaced the add with a multiply for some unknown reason. But I'm not sure why you would do that.
Anyway, how might I change your code? Fix the various things I suggested, or not. Coding style is a matter of preference, not a necessity. However, you will find that a good coding style will make your code easier to read and use in the future.
So, the only thing that is crucially important is to fix k1. I'll show you how to use a function handle to do so.
y0 = -1; % Initial Condition
% h is the increment in t
h = 0.1;
t = 0:h:100;
N = length(t);
% Pre-allocation of y_euler
y_euler = zeros(size(t));
% Initial condition gives solution at t=0.
y_euler(1) = y0;
% dy/dt
dydt = @(t,y) sqrt(t + y^2) - sqrt(t);
% Solving the equation via Euler's method
for i=1:(N-1)
k1 = dydt(t(i),y_euler(i));
y_euler(i+1) = y_euler(i) + h*k1;
end
N50=(N-1)/2;

その他の回答 (1 件)

dulanga
dulanga 2019 年 4 月 1 日
thanks for your comments John the multiplication when it was supposed to be an addition was a mistake on my part.

カテゴリ

Help Center および File ExchangeNumeric Solvers についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by