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

Question

YONGIN YI 2022-6-19

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1743325-solve-2nd-order-ode-using-matlab-ode-45-matrix

编辑： Binaya 2023-10-25

    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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

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

Answer 1

Binaya 2023-10-25

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1743325-solve-2nd-order-ode-using-matlab-ode-45-matrix#answer_1339921

编辑：Binaya 2023-10-25

在 MATLAB Online 中打开

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: