# converting python code to matlab and getting it plotted

33 visualizzazioni (ultimi 30 giorni)
Kathleen il 21 Feb 2024
Risposto: Sarthak il 26 Feb 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 CommentiMostra 1 commento meno recenteNascondi 1 commento meno recente
Kathleen il 21 Feb 2024
it was a python code. I tried to adapt it as much as possible. the goal is to get the following graphs.
Rik il 22 Feb 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?

Accedi per commentare.

### Risposta accettata

Sarthak il 26 Feb 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:
I hope the above solutions successfully resolve your query!
##### 0 CommentiMostra -2 commenti meno recentiNascondi -2 commenti meno recenti

Accedi per commentare.

### Categorie

Scopri di più su Call Python from MATLAB in Help Center e File Exchange

R2023a

### Community Treasure Hunt

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

Start Hunting!

Translated by