How do you solve a nonlinear ODE with Matlab using the finite difference approach?

Question

Yianni 2014-12-7

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/165748-how-do-you-solve-a-nonlinear-ode-with-matlab-using-the-finite-difference-approach

评论： Torsten 2014-12-8

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?

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Yianni 2014-12-7

I have attached some files to show how the coefficients were found.

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2014-12-8

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/165748-how-do-you-solve-a-nonlinear-ode-with-matlab-using-the-finite-difference-approach#answer_161681

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 个评论
显示 -1更早的评论隐藏 -1更早的评论

Torsten 2014-12-8

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.

请先登录，再进行评论。

Answer 2

Mischa Kim 2014-12-7

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/165748-how-do-you-solve-a-nonlinear-ode-with-matlab-using-the-finite-difference-approach#answer_161577

编辑：Mischa Kim 2014-12-7

在 MATLAB Online 中打开

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.

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

Yianni 2014-12-7

Thank you for this, however, unfortunately I have to use the method that I had mentioned before specifically the finite difference method. I am trying to do substitution for sin(theta) where sin(pi/6) = 0.5 (exactly) and see if it can be used in the approach I am trying to get at.

请先登录，再进行评论。