converting python code to matlab and getting it plotted

조회 수: 8 (최근 30일)
Kathleen
Kathleen 2024년 2월 21일
답변: Sarthak 2024년 2월 26일
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
Kathleen
Kathleen 2024년 2월 21일
it was a python code. I tried to adapt it as much as possible. the goal is to get the following graphs.
Rik
Rik 2024년 2월 22일
Have a read here and here. It will greatly improve your chances of getting an answer. What is the actual problem you're having?

댓글을 달려면 로그인하십시오.

채택된 답변

Sarthak
Sarthak 2024년 2월 26일
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!

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 Call Python from MATLAB에 대해 자세히 알아보기

태그

제품


릴리스

R2023a

Community Treasure Hunt

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

Start Hunting!

Translated by