solve 2nd order ODE using matlab ,ode 45, matrix,

    clc;clear;close all;

    % Define constants
 r=0.4;
 m=5;
 J=1/12*m*r^2;

 B11=0;
 B12=-3/2*m*r^2;
 B13=-1/2*m*r^2;
 B21= 3/2*m*r^2 ;
 B22= 0;
 B23= -1/2*m*r^2;
 B31= 1/2*m*r^2;
 B32= 1/2*m*r^2;
 B33= 0;

 A11=9/4*m*r^2+J;
 A12=3/2*m*r^2;
 A13=1/2*m*r^2;

 A21=3/2*m*r^2;
 A22=5/4*m*r^2+J;
 A23=1/2*m*r^2;

 A31=1/2*m*r^2;
 A32=1/2*m*r^2;
 A33=1/4*m*r^2+J;
   
    t0 = 0 ; tf = 5;   % Integration Duration
    tspan = [t0 tf];
    
    % Define initial values theta1, theta1_dot, theta2, theta2_dot
    
    theta1 = deg2rad(0); theta2 = deg2rad(0);  theta3 = deg2rad(0);  % Initial Positions (Angle)
    theta1_dot = deg2rad(0); theta2_dot = deg2rad(0); theta3_dot = deg2rad(0); % Initial Angular Velocities
    IC = [theta1,theta1_dot,theta2,theta2_dot,theta3,theta3_dot]; % Initial conditions
    
    % Define time dependent parameters (torques)
    w= pi/2;
    tau1 = @(time) 0.1*sin(w*time);  %  tau1 cos sinyali olsun
    tau2 = @(time) 0.1*sin(w*time);  %  tau2 sin sinyali olsun
    tau3 = @(time) -0.2*sin(w*time);
    
    ZERO= zeros(3,3);
    I = eye(3,3);

    x= theta1 ;
    y= theta2 ;
    z= theta3 ;
  
    x(1)= theta1_dot; 
    y(1)= theta2_dot;
    z(1)= theta3_dot;  
   
    M =[A11, A12*cos(y-x),A13*cos(z-x); A21*cos(y-x),A22,A23*cos(z-y); 
    A31*cos(z-x), A32*cos(z-y), A33];
    
    C= [B11, B12*sin(y-z)*y(1), B13*sin(z-x)*z(1); B21*sin(y-x)*x(1), B22, B23*sin(z-y)*z(1);
    B31*sin(z-x)*x(1), B32*sin(z-y)*y(1), B33];
  
    w  = pi/2; 
    
    f  =  @(time) [0.1*sin(w*time); 0.1*sin(w*time); -0.2*sin(w*time)]; 
    
    F  = M\f;	% input force vector
    
    A  = [ZERO, I ; 0 -M\C];      % state matrix
   
    Svec= [S(1); S(2); S(3); S(4); S(5); S(6)];

    SDOT= A*Svec + F;

    % Define system of first order differential equations
    
    SDOT = @(time,S) ...
     [S(4);S(5); S(6);S(7);S(8); S(9)];
    
    [time state_values] = ode45 (SDOT, tspan, IC);
    theta1 = state_values(:,1);
    theta1_dot = state_values(:,4);
    theta2 = state_values(:,2);
    theta2_dot = state_values(:,5);
    theta3 = state_values(:,3);
    theta3_dot = state_values(:,6);
    
    subplot(2,2,1)
    plot(time,rad2deg(theta1)),xlabel('time (s)'),ylabel('Joint Angle 1')   
    title(' Joint Angle 1 (Deg) vs Time')
    
    subplot(2,2,2)
    plot(time,rad2deg(theta1_dot)),xlabel('time (s)'),ylabel('Joint Velocity 1')   
    title(' Joint Velocity 1 vs Time')
    
    subplot(2,2,3)
    plot(time,rad2deg(theta2)),xlabel('time (s)'),ylabel('Joint Angle 2')   
    title(' Joint Angle 2 (Deg) vs Time')
    
    subplot(2,2,4)
    plot(time,rad2deg(theta2_dot)),xlabel('time (s)'),ylabel('Joint Velocity 2')   
    title(' Joint Velocity 2 vs Time')

    subplot(2,2,5)
    plot(time,rad2deg(theta3)),xlabel('time (s)'),ylabel('Joint Angle 3')   
    title(' Joint Angle 3 (Deg) vs Time')
    
    subplot(2,2,6)
    plot(time,rad2deg(theta3_dot)),xlabel('time (s)'),ylabel('Joint Velocity 2')   
    title(' Joint Velocity 3 vs Time')

am trying to solve the following 2nd order differential equation...[M]x''+[C]x'+[K]x=fsin(wt) where M,C and K are 3x3 matrices. I am using ODE45 to solve and ideally produce x on y-axis against time on x-axis. I am getting some errors

다음 사용 중 오류가 발생함: \

인수는 숫자형, 문자형 또는 논리형이어야 합니다.

오류 발생: untitled6 (76번 라인)

F = M\f; % input force vector

how to fix this error

or how to solve this problem?

any help is very appreciated

Answer 1

[Binaya]

Hi Yongin,

Based on your description, it seems that you want to address an error that occurs when executing the code, specifically the error message " An error occurred while using: \ The argument must be numeric, character, or logical. An error occurred: untitled6 (line 76) F = M\f; % input force vector"

Upon reviewing the provided code, it appears that the error arises at line number 76 where you try to matrix multiply inverse of “M” with “f”. The error arises because the value assigned to “f” is a function handle(defined in line number 74) and not a column vector needed as an argument for the matrix multiplication. Please consider the below suggestions for resolving the error:

1. Pass a value as an argument while calling “f”:

F  = M\f(10);	% input force vector 

2. Change “f” from a function handle to a constant column vector.

3. Please assign values to the "time" variable, which is subsequently used to define the variable "f".

Please note that the rest of the code appears to be correct and does not contribute to the aforementioned issue.

I hope this helps.

Regards

Binaya