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
Best Answer