I have attached some files to show how the coefficients were found.
How do you solve a nonlinear ODE with Matlab using the finite difference approach?
4 views (last 30 days)
Show older comments
I have done this with a linear ODE which had the equation x''+(c/m)*x'+(g/L)*x = 0 where x(0) = pi/6 and x'(0) = 0
%Method 2: Numerical Solution Using the Finite Difference Approach
clear all,close all
m = 1; %Mass of Pendulum (kg)
c = 2; %Friction Coefficient (kg/s)
L = 1; %Length of Pendulum Arm (m)
g = 10; %Gravitational Acceleration (m/s^2)
Nt = 101; %Step Size of time
ti = 0; %Initial time (sec)
tf = 10; %Final time (sec)
t = linspace(ti,tf,Nt); %Time vector (sec)
x1 = pi/6; %Initial Position (radians)
v1 = 0; %Initial Velocity (radians/s)
N = Nt-1; dt = (tf-ti)/N; %dt is the change of t over N which is the step size
%Evaluated Equation Coefficients with Starting Points
% xn is the Angular Position (degrees) of Case 1, Method 2
a = 1 + c*dt/(2*m);
b = 2 - g*dt*dt/L;
d = c*dt/(2*m) - 1;
xn = zeros(1,Nt); xn(1) = x1; xn(2) = xn(1) + v1*dt;
%Loop Over Remaining Discrete Time Points
for i = 2:N
xn(i+1) = b*xn(i)/a + d*xn(i-1)/a;
end
figure(1)
plot(t,xn*180/pi),grid on
title('Linear Model Behavior for Case 1')
xlabel('Time (sec)'), ylabel('Angular Position (degrees)')
This file represents a solution using a finite difference approach for a linear ODE. I would like to do the same with a nonlinear ODE specifically x''+(c/m)*x'+(g/L)*sin(x) = 0 where x(0) = pi/6 and x'(0) = 0. (THE DIFFERENCE IS THE USE OF THE SIN FUNCTION). How can this be done?
Accepted Answer
Torsten
on 8 Dec 2014
All you have to do is replace
b = 2 - g*dt*dt/L;
by
b=2;
and write the new recurrence as
xn(i+1) = b*xn(i)/a + d*xn(i-1)/a - g/L*sin(xn(i))/a.
Best wishes
Torsten.
1 Comment
Torsten
on 8 Dec 2014
Sorry, should read
xn(i+1) = b*xn(i)/a + d*xn(i-1)/a - g/L*dt^2*sin(xn(i))/a.
Best wishes
Torsten.
More Answers (1)
Mischa Kim
on 7 Dec 2014
Edited: Mischa Kim
on 7 Dec 2014
Yianni, unless you want/have to re-invent the wheel use one of MATLAB's integrators, e.g., ode45
function test_ode()
m = 1;
c = 2;
L = 1;
g = 10;
param = [c; m; g; L];
IC = [pi/6; 0];
tspan = linspace(0,10,100);
[T,X] = ode45(@my_de,tspan,IC,[],param);
plot(T,X(:,1),T,X(:,2))
end
function dx = my_de(t,x,param)
c = param(1);
m = param(2);
g = param(3);
L = param(4);
r = x(1);
v = x(2);
dx = [v;...
-(c/m)*v - (g/L)*sin(r)];
end
Put both functions in one and the same file and save it as test_ode.m.
See Also
Categories
Find more on Assembly in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!