euler method first order

25 vues (au cours des 30 derniers jours)
dulanga
dulanga le 1 Avr 2019
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

Réponse acceptée

John D'Errico
John D'Errico le 1 Avr 2019
Modifié(e) : John D'Errico le 1 Avr 2019
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;

Plus de réponses (1)

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

Community Treasure Hunt

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

Start Hunting!

Translated by