Euler Method Solving by Differential Equation - Lecture Notes | BME 201, Study notes of Biology

Outline Feb 28 Material Type: Notes; Professor: Sawicki; Class: Computer Methods in Biomedical Engineering; Subject: Biomedical Engineering; University: North Carolina State University; Term: Spring 2011;

Typology: Study notes

2010/2011

Uploaded on 04/06/2011

goalie4eva
goalie4eva 🇺🇸

31 documents

1 / 5

Toggle sidebar

This page cannot be seen from the preview

Don't miss anything!

bg1
Lecture 13: Monday February 28, 2011:
Announcements:
Nothing is due today .
Case 2 resubmission is due Wednesday, March 2. Be sure to submit you’re your
original and the new version of your solution. Please take advantage of Prof.
Sawicki’s Office Hours on M’s and T’s from 4-5PM 4212C EB3. Ben’s hours are
still being held Thursday afternoons.
SH6 and Case 3 will be posted by the end of the week. These will be due Wed.
March 16th (SH6) and Monday March 21st (Case 3).
SL5 is not due until after Spring Break (Thursday March 17th).
Lab this Thursday will be a demonstration of a powered ankle exoskeleton to
give you a flavor of what is to come following spring break.
I am still grading the exam and hope to return in Monday March 14th. Thanks for
your patience.
Reading for Next Time:
Reading on Functions and Program Development: Palm Chapter 3.1-3.3, Palm Chapter
4.1 Reading on Debugging Palm Chapter 4.8
*Today’s Goals:
Review Euler’s Method (e.g. similar to SL5). Begin discussing functions- MATLAB built-
in and user-defined.
*I REMEMBER: turn on your diary to save your Command Window
>>diary LectureFeb28.m
*II. Euler’s Method for solving differential equations, revisited.
There are many differential equations that cannot be solved analytically (i.e. on pencil
and paper). For these we need a computer to march through time-
Explain Euler’s method on whiteboard.
As an example, derive simple pendulum equations of motion and make them into two
first order equations.
Step 1. deltaY/deltaT=somediffEquation
Step 2. y1 = y0 + deltay
Step 3. deltay =somediffEquation*dt
%Euler's Method Example:
%Integrating equations of motion describing pendulum mechanics.
%PROBLEM DESCRIPTION: Start pendulum at angle theta=pi/2 radians
%and angular velocity omega=0 rad/s. Find the time when the
%pendulum passes through 0 radians. Your code should be flexible
%enough to handle pendula of different lengths and masses as well
%as different initial conditions.
clear all
pf3
pf4
pf5

Partial preview of the text

Download Euler Method Solving by Differential Equation - Lecture Notes | BME 201 and more Study notes Biology in PDF only on Docsity!

Lecture 13: Monday February 28, 2011:

Announcements:

 Nothing is due today .

 Case 2 resubmission is due Wednesday, March 2. Be sure to submit you’re your

original and the new version of your solution. Please take advantage of Prof.

Sawicki’s Office Hours on M’s and T’s from 4-5PM 4212C EB3. Ben’s hours are

still being held Thursday afternoons.

 SH6 and Case 3 will be posted by the end of the week. These will be due Wed.

March 16th^ (SH6) and Monday March 21st^ (Case 3).

 SL5 is not due until after Spring Break (Thursday March 17th).

 Lab this Thursday will be a demonstration of a powered ankle exoskeleton to

give you a flavor of what is to come following spring break.

 I am still grading the exam and hope to return in Monday March 14th. Thanks for

your patience.

Reading for Next Time:

Reading on Functions and Program Development: Palm Chapter 3.1-3.3, Palm Chapter

4.1 Reading on Debugging Palm Chapter 4.

*Today’s Goals:

Review Euler’s Method (e.g. similar to SL5). Begin discussing functions- MATLAB built-

in and user-defined.

*I REMEMBER: turn on your diary to save your Command Window

>>diary LectureFeb28.m

*II. Euler’s Method for solving differential equations, revisited.

There are many differential equations that cannot be solved analytically (i.e. on pencil

and paper). For these we need a computer to march through time-

Explain Euler’s method on whiteboard.

As an example, derive simple pendulum equations of motion and make them into two

first order equations.

Step 1. deltaY/deltaT=somediffEquation

Step 2. y1 = y0 + deltay

Step 3. deltay =somediffEquation*dt

%Euler's Method Example:

%Integrating equations of motion describing pendulum mechanics.

%PROBLEM DESCRIPTION: Start pendulum at angle theta=pi/2 radians

%and angular velocity omega=0 rad/s. Find the time when the

%pendulum passes through 0 radians. Your code should be flexible

%enough to handle pendula of different lengths and masses as well

%as different initial conditions.

clear all

%Set pendulum and environment parameters

gravity=9.81; %m/s^

length=1; %m pendulum length

mass=10; %kg pendulum mass

%Set dt and initial conditions for all variables

dt=.0005; %s %Set time step size for numerical integration.

theta0=pi/2; %rad -Positve is counterclockwise from bottom

omega0=0; %rad/s -Starts at rest

%Set event condition to control when simulation should stop.

%Use this in while loop control.

thetalimit=-pi/2; %rad

time(1)=0; %s

theta(1)=theta0;%rad

omega(1)=omega0; %rad/s

index=1; %initialize index for time, theta and omega array

while theta(index) > thetalimit %Check if theta is at limit yet

%1.Update the value of pendulum angular velocity based on

%Newtons' Laws

omega(index+1)= omega(index)-gravitysin(theta(index))/lengthdt;

%2.Update the value of pendulum angle based on

theta(index+1)= theta(index)+ omega(index)*dt;

%3.Advance time one timestep

time(index+1)=time(index)+dt;

%4.Increase array index by 1 to move to next location for storage

index=index+1;

end

%5.Print the time when condition is satisfied

tfinal=time(end)

thetadeg=theta*180/pi; %convert radians to degrees

omegadeg=omega*180/pi; %convert radians to degrees

%6.Create subplot for angle and angular velocity

figure(1)

subplot(2,1,1)

plot(time,thetadeg)

title('Pendulum Angle and Angular Velocity vs. Time')

ylabel('Angle (deg)')

grid on

subplot(2,1,2)

plot(time,omegadeg)

xlabel('Time (sec)')

ylabel('Angular Velocity (deg/s)')

grid on

Here is how I could use the function computeThis1 - by ‘calling’ it from a separate script

file (e.g. main program.m).

%define three variables a=4, b=3, c=7; %call function that takes 3 variables and returns 3 variables. [d, e, f] = computeThis1(a, b, c);

The function does not change the values of a, b, or c but instead, it returns 3 values that

I chose to save as d, e, and f. I could have chosen to overwrite a, b, c, but in this case I

want to save the original values.

NOTE WELL: The variable names in the function definition file do not have to match the

variable names that I use when ‘calling’ the function. The order that the variables are

passed in will be maintained, that is, x=a, y=b, and z=c inside the function. And on the

return, d=x, e=y, and f=z.

It is also possible to nest functions (i.e. subfunctions) within other functions

(primary/parent functions). For example:

function [x, y, z] = computeThis2(x, y, z) % a simple function that takes 3 input values, % and returns the same 3 values that may or may not have been updated in %the function body. This is the parent/primary function. if x==y %use == to compare for equivalence y = z x = x^2; elseif x < y %less than y = x+1; else %if previous statements are false... [x, y, z]= subElse(x, y, z) end % terminate code block end % terminate parent function %A subfunction that can only be used by functions in this file function [l, m, n] = subElse(l, m, n) m = l^2; l = n; end

In the above example, when computeThis2 uses the subElse function, and we overwrite

the original values for x, y, and z since these are being returned inside the parent

function.

SOME RULES FOR FUNCTIONS:

*Rule 1. To build a user-defined function you must have consistency between the m-file

name and the function name. Don’t use a name already in use in MATLAB. Use the

exist command to check if you are unsure.

>> help exist

EXIST Check if variables or functions are defined.

0 if A does not exist

1 if A is a variable in the workspace

2 if A is an M-file on MATLAB's search path.

*Rule 2. You must follow proper syntax in the function definition line. That is, the function

output arguments (what it computes) go on the left of the = in [ ] brackets and the

function input arguments (what it works on) go on the right of the = in () parentheses.

For example the basic structure of the user-defined function FunctionName.m is as

follows:

The first line of the file FunctionName.m must contain-

function[output_args] = FunctionName(input_args)

then,

%FunctionName.m Summary of this function goes here

%Detailed explanation of function goes here

*CODE GOES HERE!!!!!

end

*Rule 3. Variables inside the function file are called local variables and cannot be

accessed except from within the function itself (i.e. they are protected).

*Rule 4. When calling a function, only the order of the inputs/arguments ( ) and

outputs [ ] matter, not the names.

*Rule 5. Functions can have any number of inputs or outputs (including no inputs or

outputs). For example, sometimes a function just returns a plot or screen display and

has no outputs. Or sometimes a function just initializes some local variables and has no

inputs. See pages 121 and 122 in Palm for other examples.

*Rule 6. NEVER define global variables! i.e. Ignore Palm top of pg. 124.

*Group Exercise – Q2: How could we modify your pendulum simulation script file so

that it is a user-defined function called PendSim. It should take

length,grav,tfinal,dt,ang0,angvel0 as inputs and return time,ang,angvel.

See posted code on Moodle for solution.

*NEXT TIME: More on functions! Including function handles, function functions, calling

functions, anonymous functions etc.