Effacer les filtres
Effacer les filtres

converting python code to matlab and getting it plotted

26 vues (au cours des 30 derniers jours)
Kathleen
Kathleen le 21 Fév 2024
Réponse apportée : Sarthak le 26 Fév 2024
I been trying to get this python code to work in matlab. Besides converting issues I probably have a lot of other mistakes in it as I am just a beginner.
% constants
L = 100; % length of bar in km
dx = 1; % spatial grid spacing km
beta = 4; % shear wave velocity in km/s
dt = 0.1; % time spacing in seconds
T = 5; % total time in sec
N = int(T/dt); % number of time steps
tmax = 30; % wave run time sec
t = [0:dt:tmax]; %time vec
nt = length (t);
% Define initial conditions
u = zeros(int(L/dx)+1); % displacement array
u(int(50e3/dx)) = 1; % source at u(50 km)
% intialize 3 displacement vectors
%u3 = new solution
%u2 = previous solution
%u1 = previous previous solution
% displacement arrays
u1 = zeroes (1, 101);
u2 = u1;
u3 = u1;
% boundary conditions
u(0) = 0;
u(-1) = 0;
% Define finite difference coefficients
A = (beta*dt/dx)^2
B = 2 - 2*A
% Solve wave equation using finite differences
for n = range(1,N)
u(1:-1) = B*u(1:-1) - u(1:-1) + A*(u (2: + u :-2))
u(int(50e3/dx)) = np.sin(2*np.pi*n*dt/5)
end
% Plot solution
figure(1); hold on; grid on; axis equal
xline(0,'-','linewidth',2); yline(0,'-','linewidth',2) %adding x and y line
set(0,'DefaultLineLineWidth',5,'DefaultAxesFontSize',12)
x = linspace(0, L, int(L/dx)+1)
plot(x, u)
xlabel('Distance (m)')
ylabel('Displacement')
show()
% initial condition
u3(50) = sin(pi*tlen/5)^2
% boundary conditions
u3(0) = u3(1) % stress-free boundary
u3(L) = 0 % fixed boundary
% Define plotting function
figure plot_u(u, t):
x = arange(0, L+dx, dx)
plot(x, u)
title('t = {t:.2f} s')
xlabel('x (km)')
ylabel('u (m)')
ylim(-0.1, 0.1)
show()
%%% Solve PDE using finite differences
t = 0
if t <= tlen
t = dt
for i [1, L]
rhs = beta^2 * (u2[i+1] - 2*u2[i] + u2[i-1]) / dx^2
u3[i] = 2*u2[i] - u1[i] + dt^2 * rhs
% Apply boundary conditions
u3(0) = u3(1) % stress-free boundary
u3(L) = 0 % fixed boundary
if t <= tlen
u3(50) = sin(np.pi*t/tlen)^2
% Update displacement arrays
u1 = copy(u2)
u2 = copy(u3)
% Plot every 4 s
if t % 4 < dt:
plot_u(u2, t)
end
end
end
end
%% Calculate velocities of pulses
t1 = 1.5 % time to measure velocity of first pulse
t2 = 12.5 % time to measure velocity of second pulse
x1 = 25 % distance to measure velocity of first pulse
x2 = 75 % distance to measure velocity of second pulse
vel1 = (u2(x1+1,t1) - u2(x1-1,t1)) / (2*dx*dt)
vel2 = (u2(x2+1,t2) - u2(x2-1,t2)) / (2*dx*dt)
fprintf('Velocity of first pulse: {vel1:.2f} km/s')
fprintf('Velocity of second pulse: {vel2:.2f} km/s')
figure(2) % displacement at endpoints as a function of time
x_endpoints = [0, L]
u_endpoints = u2*[0,':'], u2*[L,':']
for i % in range(2):
plot(arange(0, tlen+dt, dt), u_endpoints[i])
xlabel('t (s)')
ylabel('f u(x_endpoints[i] km) (m)')
show()
% Plot displacement at crossing
(initialize t, dx, dt, tlen, beta and u1, u2, u3 arrays)
while t <= 33
t = t + dt;
for i = 2:99
rhs = beta^2 * (u2(i+1) - 2*u2(i) + u2(i-1)) / dx^2;
u3(i) = dt^2 * rhs + 2*u2(i) - u1(i);
end
u3(1) = u3(2); % stress-free boundary
u3(100) = 0; % fixed boundary
if t <= tlen
u3(50) = sin(pi * t/tlen)^2;
end
u1 = u2;
u2 = u3;
% output u2 at desired times
end
  3 commentaires
Kathleen
Kathleen le 21 Fév 2024
it was a python code. I tried to adapt it as much as possible. the goal is to get the following graphs.
Rik
Rik le 22 Fév 2024
Have a read here and here. It will greatly improve your chances of getting an answer. What is the actual problem you're having?

Connectez-vous pour commenter.

Réponse acceptée

Sarthak
Sarthak le 26 Fév 2024
Hi Kathleen,
As per my understanding, there is no direct way to convert python code to MATLAB.
However you can try the following approaches:
Please refer to the following MATLAB answers post with similar issue: https://www.mathworks.com/matlabcentral/answers/416129-how-to-convert-python-code-into-matlab
I hope the above solutions successfully resolve your query!

Plus de réponses (0)

Catégories

En savoir plus sur Call Python from MATLAB dans Help Center et File Exchange

Produits


Version

R2023a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by