Solving a differential equation using ode45

B1, B2, and B3 are 4-D double, and intial condition are phi = 0 at s = 0, d(phi)/ds = 0 at s = 0. I am trying to solve this equation using ode45. EQUATION IS B1*phi + B2*(d(phi)/ds) - B3*(d^2(phi)/ds^2) = 0.

my code

% Define the linear coefficient matrices B at the specific s value
for i = 1:n_node
    B(:,:,:,i) = coefficient_linear(nmp_order, GPr, GPt, GPWr, GPWt, NGPr, NGPt, s(i));
end
% Initialize variables
B1 = B(1,:,:,:);    % Coefficient B1
B2 = B(2,:,:,:);    % Coefficient B2
B3 = B(3,:,:,:);    % Coefficient B3
% Define the anonymous function for the equation
equation = @(s, phi_dphi) [phi_dphi(2); -B1.*phi_dphi(1) + B2.*phi_dphi(2) - B3.*(phi_dphi(2).^2)];
% Define the initial conditions
phi0 = 0;              % Initial condition for phi
dphi_ds0 = 0;          % Initial condition for d(phi)/ds
% Set the options for ode45 (optional)
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-6);
% Solve the equation using ode45
[t, phi_dphi] = ode45(equation, [S_0, S_L], [phi0; dphi_ds0], options); % error at this line of code 
% Extract the solution for phi and d(phi)/ds
phi = phi_dphi(:, 1);
dphi_ds = phi_dphi(:, 2);
% Display the results
disp('Solution:')
disp('---------')
disp('s       phi       d(phi)/ds')
disp([t, phi, dphi_ds]);

I have bold the line of code and the error is

Error using vertcat

Dimensions of arrays being concatenated are not consistent.

[KSSV]

Check the dimensions of :

[S_0, S_L]
[phi0; dphi_ds0]

[GOBIND KUMAR]

[S_0, S_L] = [0,6] % dimension of the tube

[phi0; dphi_ds0] = [0;0] % initial condition

does the dimension of these 2 needs to be 4x1 column vector?

Answer 1

[Sam Chak]

Hi @GOBIND KUMAR,

The code for the ode45 part can run. No issue with that.

So, I think the problem lies in the computation of the . Try fixing that.

S_0 = 0;
S_L = 20;

% Initialize variables
B1 =  1;    % Coefficient B1
B2 = -2;    % Coefficient B2
B3 =  0;    % Coefficient B3

% Define the anonymous function for the equation
equation = @(s, phi_dphi) [phi_dphi(2); -B1.*phi_dphi(1) + B2.*phi_dphi(2) - B3.*(phi_dphi(2).^2)];

% Define the initial conditions
phi0     = 1;       % Initial condition for phi
dphi_ds0 = 0;       % Initial condition for d(phi)/ds

% Set the options for ode45 (optional)
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-6);

% Solve the equation using ode45
[t, phi_dphi] = ode45(equation, [S_0, S_L], [phi0; dphi_ds0], options); % error at this line of code 

% Extract the solution for phi and d(phi)/ds
phi = phi_dphi(:, 1);
dphi_ds = phi_dphi(:, 2);

% Test plot
plot(t, phi), grid on

% Display the results
disp('Solution:')

Solution:

disp('---------')

---------

disp('s       phi       d(phi)/ds')

s       phi       d(phi)/ds

disp([t, phi, dphi_ds]);

         0    1.0000         0
    0.0126    0.9999   -0.0125
    0.0252    0.9997   -0.0246
    0.0379    0.9993   -0.0365
    0.0505    0.9988   -0.0480
    0.0883    0.9963   -0.0808
    0.1261    0.9927   -0.1112
    0.1639    0.9879   -0.1391
    0.2017    0.9822   -0.1649
    0.2396    0.9755   -0.1886
    0.2775    0.9679   -0.2103
    0.3154    0.9596   -0.2301
    0.3533    0.9505   -0.2482
    0.3927    0.9404   -0.2651
    0.4320    0.9297   -0.2804
    0.4713    0.9184   -0.2942
    0.5106    0.9066   -0.3064
    0.5513    0.8938   -0.3177
    0.5921    0.8807   -0.3275
    0.6329    0.8672   -0.3361
    0.6736    0.8533   -0.3435
    0.7159    0.8386   -0.3499
    0.7583    0.8237   -0.3552
    0.8006    0.8086   -0.3595
    0.8429    0.7933   -0.3628
    0.8869    0.7773   -0.3653
    0.9309    0.7611   -0.3670
    0.9750    0.7450   -0.3678
    1.0190    0.7288   -0.3678
    1.0648    0.7119   -0.3671
    1.1107    0.6951   -0.3658
    1.1566    0.6784   -0.3638
    1.2024    0.6617   -0.3613
    1.2503    0.6445   -0.3581
    1.2982    0.6275   -0.3544
    1.3461    0.6106   -0.3503
    1.3939    0.5939   -0.3458
    1.4440    0.5767   -0.3408
    1.4941    0.5598   -0.3353
    1.5442    0.5431   -0.3297
    1.5943    0.5268   -0.3237
    1.6468    0.5099   -0.3173
    1.6994    0.4934   -0.3106
    1.7519    0.4773   -0.3039
    1.8044    0.4615   -0.2970
    1.8597    0.4453   -0.2896
    1.9149    0.4295   -0.2822
    1.9702    0.4141   -0.2747
    2.0254    0.3992   -0.2672
    2.0837    0.3838   -0.2594
    2.1419    0.3689   -0.2515
    2.2002    0.3545   -0.2437
    2.2585    0.3405   -0.2360
    2.3202    0.3262   -0.2280
    2.3818    0.3124   -0.2200
    2.4435    0.2991   -0.2122
    2.5052    0.2862   -0.2046
    2.5708    0.2731   -0.1966
    2.6363    0.2604   -0.1888
    2.7019    0.2483   -0.1812
    2.7674    0.2367   -0.1739
    2.8375    0.2248   -0.1662
    2.9075    0.2134   -0.1588
    2.9775    0.2025   -0.1516
    3.0475    0.1922   -0.1447
    3.1228    0.1815   -0.1375
    3.1980    0.1715   -0.1306
    3.2733    0.1619   -0.1240
    3.3485    0.1528   -0.1177
    3.4299    0.1435   -0.1111
    3.5114    0.1347   -0.1048
    3.5928    0.1264   -0.0989
    3.6742    0.1186   -0.0932
    3.7632    0.1105   -0.0873
    3.8522    0.1030   -0.0818
    3.9412    0.0960   -0.0766
    4.0302    0.0894   -0.0716
    4.1287    0.0826   -0.0665
    4.2273    0.0763   -0.0617
    4.3258    0.0704   -0.0572
    4.4243    0.0650   -0.0530
    4.5356    0.0593   -0.0486
    4.6469    0.0542   -0.0446
    4.7581    0.0494   -0.0408
    4.8694    0.0451   -0.0374
    4.9996    0.0404   -0.0337
    5.1298    0.0363   -0.0304
    5.2599    0.0325   -0.0273
    5.3901    0.0291   -0.0246
    5.5343    0.0258   -0.0219
    5.6785    0.0228   -0.0194
    5.8228    0.0202   -0.0172
    5.9670    0.0178   -0.0153
    6.1096    0.0158   -0.0136
    6.2523    0.0140   -0.0120
    6.3949    0.0123   -0.0107
    6.5375    0.0109   -0.0095
    6.6854    0.0096   -0.0084
    6.8332    0.0084   -0.0074
    6.9811    0.0074   -0.0065
    7.1289    0.0065   -0.0057
    7.2854    0.0057   -0.0050
    7.4419    0.0049   -0.0044
    7.5984    0.0043   -0.0038
    7.7549    0.0038   -0.0033
    7.9229    0.0032   -0.0029
    8.0909    0.0028   -0.0025
    8.2590    0.0024   -0.0021
    8.4270    0.0021   -0.0018
    8.6097    0.0018   -0.0016
    8.7924    0.0015   -0.0013
    8.9750    0.0013   -0.0011
    9.1577    0.0011   -0.0010
    9.3588    0.0009   -0.0008
    9.5599    0.0007   -0.0007
    9.7610    0.0006   -0.0006
    9.9621    0.0005   -0.0005
   10.1865    0.0004   -0.0004
   10.4108    0.0003   -0.0003
   10.6352    0.0003   -0.0003
   10.8595    0.0002   -0.0002
   11.1137    0.0002   -0.0002
   11.3679    0.0001   -0.0001
   11.6220    0.0001   -0.0001
   11.8762    0.0001   -0.0001
   12.1693    0.0001   -0.0001
   12.4625    0.0001   -0.0000
   12.7557    0.0000   -0.0000
   13.0488    0.0000   -0.0000
   13.3945    0.0000   -0.0000
   13.7401    0.0000   -0.0000
   14.0858    0.0000   -0.0000
   14.4315    0.0000   -0.0000
   14.8505    0.0000   -0.0000
   15.2696    0.0000   -0.0000
   15.6887    0.0000   -0.0000
   16.1077    0.0000   -0.0000
   16.6077    0.0000   -0.0000
   17.1077    0.0000   -0.0000
   17.6077    0.0000   -0.0000
   18.1077    0.0000   -0.0000
   18.5808    0.0000   -0.0000
   19.0539    0.0000   -0.0000
   19.5269    0.0000   -0.0000
   20.0000    0.0000   -0.0000

[GOBIND KUMAR]

sir, here you have taken B1,B2, B3 as 1x1 scalar but it is a 4D Double matrix. while S_0,S_L,phi,dphi_ds are 1x1 scalar

[Sam Chak]

Hi @GOBIND KUMAR

Thanks for your clarification. How do you interpret this ODE?

equation = @(s, phi_dphi) [phi_dphi(2); -B1.*phi_dphi(1) + B2.*phi_dphi(2) - B3.*(phi_dphi(2).^2)];

Does it look this?

How exactly do you perform the 4-D matrix multiplication? If you can show us the math, it will help to resolve the issue.

[GOBIND KUMAR]

sir, I am not sure about the multiplication part using B1,B2,B3.

[Sam Chak]

Hi @GOBIND KUMAR

I'm not too familiar with some of your math notations. But I think in vector field form, the ODEs should be expressed as follows:

If that is true, then your original ode45 code above was incorrect.

If you are unsure about the matrix multiplication B1, B2, B3, perhaps you can plot out B1, B2, B3. You need to at least provide info about the B's. The ode45 part does not look like a problem yet.