Solve first-order ordinary differential equations

Hello,

I need to solve the following system of differential equations:

where V(x,y) is a matrix that does not necessarily have a determined value for the integration points (I believe that it is necessary to interpolate to generate such points).

I implemented the following code, but I do not trust the results. I expected a greater variation in the phi angle. Could someone please help validate my code?

function [time, x, y, phi] = my_solution

    load('Data.mat') % attached file
    
    x = linspace(0,50e-3,500);
    y = linspace(0,50e-3,500);
    [Xm, Ym] = meshgrid(x,y); % grid
    
    % Vm is the matrix of velocitys
    grad_xm = gradient(Vm,1); % partial derivative with respect to x
    grad_ym = gradient(Vm,2); % partial derivative with respect to y
    
    % [S] = [s(1) s(2) s(3)]' = [x y phi]'
    % sdot = [S'] = [x' y' phi']' = [V(x,y)cos(phi)   V(x,y)sin(phi)   sin(phi)*grad_xm(x,y)-cos(phi)*grad_ym(x,y)]'

    % functions for interpolation at the integration points.
    V = @(t, s) interp2(Xm,Ym,Vm,s(1),s(2));
    grad_x = @(t, s) interp2(Xm,Ym,grad_xm,s(1),s(2));
    grad_y = @(t, s) interp2(Xm,Ym,grad_ym,s(1),s(2));
    
    sdot = @(t,s)[V(t,s)*cos(s(3)); ...
        V(t,s)*sin(s(3)); ...
        sin(s(3))*grad_x(t,s) - cos(s(3))*grad_y(t,s)];
    
    tspan = [0,100e-6];
    IC = [20e-3 0 pi/3]
    
    [time, state_values] = ode45(sdot,tspan,IC);
    
    x = state_values(:,1);
    y = state_values(:,2);
    phi = state_values(:,3);
    
    figure
    subplot(311)
    plot(time,x)
    ylabel('x')
    subplot(312)
    plot(time,y)
    ylabel('y')
    subplot(313)
    plot(time,phi*180/pi)
    ylabel('\phi')
    xlabel('Time')

end

Answer 1

[Bjorn Gustavsson]

Your calculation of grad_xm and grad_ym looks peculiar to me. Since you seem to use x and y for the coordinates of Vm, it seems like you'd want to calculate the gradients like this:

[grad_xm,grad_ym] = gradient(Vm,x,y);

If I use that then I get values of phi that varies between 60 and 45 degrees.

Your calculation of grad_xm and grad_ym uses the spacing of 1 and 2 respectively for the spatial coordinates when you intend to have 50e3/499 as spacing.

HTH

[Juliano Scholz Slongo]

Bjorn Gustavsson,

You're absolutely right. I was using the function in a wrong way.

Your help was effective.

Thank you so much.

[Bjorn Gustavsson]

My pleasure.

(These things are true: To err is human, to realy screw up you need a computer. Today I found an error in my code that's too embarrasing to share - just know that you (sadly) most likely will continue to make bugs - we just have to keep on fighting "them"...)