58:110 Computer-Aided Engineering Spring 2005

MATLAB EXAMPLES

Initial-value problems Department of Mechanical and industrial engineering January 2005

Some useful functions

ode23 ode23s ode23t ode23tb ode45 ode113 ode15s

An explicit one-step Runge-Kutta low order(2 nd -3 rd )solver An implicit one-step modified Rosenbrock solver of order two An implicit one-step trapezoidal rule using a free interpolate An implicit trapezoidal rule followed by a backward differentiation rule of order two An explicit one-step Runge Kutta medium-order (4 th -5 th ) solver A multistep Adams-Bashforth-Moulton PECE solver of varying order (1 st –13 th ) An implicit, multistep numerical differentiation solver of varying order (1 st -5 th )

IVP

The initial value problems for 1 st order differential equation are defined by the following standard form:

x

     

is given

Example:

x x

1

x

(0)



0 The exact solution is:

x e

1

x

In the following slides, we will try to solve this problem by Euler and Runge-Kutta methods in MATLAB. The ODE function was defined as M-file (ex1_odefun.m): % IVP function for EX1 function xdot=ex1_odefun(t,x) xdot=x+1;

IVP –Euler method

Use Euler method to solve the IVP example, step size sets to 0.01.

Write the M-file as follows (ivp-euler.ex1.m): % IVP solution - Euler method % example 1: x'=x+1, x(0)=0, solve x(1) % the range of t [a,b] a=0.0; b=1.0; % the initial value of x x=0.0; %the number of steps n=100; %set the matrix x0=zeros(n,1); % set the step size h=(b-a)/100; %set the initial value of t t=a; for k=1:n % Record the values of x in each step x=x+h*ex1_odefun(t,x); x0(k)=x; t=t+h; end %the final result of x x0(100) Run this M-file: >> ivp_euler_ex1 ans = 1.7048

IVP –Runge-Kutta method

To use Runge-Kutta method, here we only try ode23 and ode45 (step size =0.1): >> tspan=linspace(0,1,11); >> x0=[0;0]; >> [t,x]=ode23(@ex1_odefun,tspan,x0); >> x(11) ans = 1.7182

>> [t,x]=ode45(@ex1_odefun,tspan,x0); >> x(11) ans = 1.7183

The results show the higher order solver ode45 got more accurate result than ode23 and Euler method.

IVP - Plot

To show the accuracy of different methods, we need plot the solution together, the following M file is for plotting the results of exact, Euler and 4 th Runge-Kutta methods (ivp-plot_ex1.m): % IVP solution - plot % example 1 x'=x+1 % Record the values of x in each step % the range of t [a,b] a=0.0; x1(k)=x; x=x+h*ex1_odefun(t,x); b=1.0; % the initial value of x x=0.0; %the number of steps n=10; %set the matrix x1=zeros(n,1); %the step size h=(b-a)/10; t=t+h; end %plot the results t=linspace(0,1,11); %Runge-Kutta method rx0=[0;0]; [t,x2]=ode45(@ex1_odefun,t,rx0); %exact solution x0=exp(t)-1; %set the initial value of t t=a; plot(t,x0,t,x1,'.',t,x2,'o') xlabel('t') %Euler method - Iterations for k=1:n+1 ylabel('x')

IVP - Plot

The following figures show the results by different step size h. The solid line is the exact solution, the solid dots represent the result of Euler method, the circles represents the results of 4 th order Runge-Kutta method.

h = 0.1

h = 0.01

IVP-Summary

Compare the results of different methods (all keep 4 digits after the radix point) in the following table (x(1)):

Method

Exact Euler (h = 0.1) 4 th Runge-Kutta (h = 0.1) Euler (h = 0.01) 4 th Runge-Kutta (h = 0.01)

Solution

1.7183

1.5937

1.7183

1.7084

1.7183