Using numerical methods to plot solution to first-order nonlinear differential equation in Matlab-CodePudding

I have a question about plotting x(t), the solution to the following differential equation knowing that dx/dt equals the expression below. The value of x is 0 at t = 0.

syms x

dxdt = -(1.0*(6.84e 45*x^2   5.24e 32*x - 2.49e 42))/(2.47e 39*x   7.12e 37)

I want to plot the solution of this first-order nonlinear differential equation. The analytical solution involves complex numbers so that's not relevant because this equation models a real-life process, but Matlab can solve the equation using numerical methods and plot it. Can someone please suggest how to do this?

CodePudding user response：

in matlab try this

tspan = [0 10];
x0 = 0;
[t,x] = ode45(@(t,x) -(1.0*(6.84e 45*x^2   5.24e 32*x - 2.49e 42))/(2.47e 39*x   7.12e 37), tspan, x0);
plot(t,x,'b')

i try it and i got this

hope that help you.

CodePudding user response：

I have written an example for how to use Python with SymPy and matplotlib. SymPy can be used to calculate both definite and indefinite integrals. By calculating the indefinite integral and adding a constant to set it to evaluate to 0 at t = 0. Now you have the integral, so just a matter of plotting. I would define an array from a starting point to an endpoint with 1000 points between (could likely be less). You can then calculate the value of the integral with the constant at each time point, which can then be plotted with matplotlib. There are plenty of other questions on how to customize plots with matplotlib.

This displays a basic plot of the indefinite integral of the function dxdt with assumption of x(t) = 0. Variation of the tuple when running Plotting() will set what range of x values to plot. This is set to plot 1000 data points between the minimum and maximum values set when calling the function.

For more information on customizing the plot, I recommend matplotlib documentation. Documentation on the integral can be found in SymPy documentation.

import pandas as pd
from sympy import *
from sympy.abc import x
import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np

def Plotting(xValues, dxdt):
    # Calculate integral
    xt = integrate(dxdt,x)
    # Convert to function
    f = lambdify(x, xt)
    C = -f(0)

    # Define x values, last number in linspace corresponding to number of points to plot
    xValues = np.linspace(xValues[0],xValues[1],500)
    yValues = [f(x) C for x in xValues]
    
    # Initialize figure
    fig = plt.figure(figsize = (4,3))
    ax = fig.add_axes([0, 0, 1, 1])
    
    # Plot Data    
    ax.plot(xValues, yValues)
    plt.show()
    plt.close("all")

# Define Function
dxdt = -(1.0*(6.84e45*x**2   5.24e32*x - 2.49e42))/(2.47e39*x   7.12e37)

# Run Plotting function, with left and right most points defined as tuple, and function as second argument
Plotting((-0.025, 0.05),dxdt)