-1

I have an equation like this:

dy/dx = a(x)*y + b

where a(x) is a non-constant (a=1/x) and b is a vector (10000 rows).

How can I solve this equation?

Sardar Usama
  • 19,536
  • 9
  • 36
  • 58
Ari
  • 15
  • 4
  • 2
    You don't need MATLAB to solve this: you can use an [integrating factor](https://en.wikipedia.org/wiki/Integrating_factor) to find a solution `y(x) = b*x*ln(x) + c*x` for some constant c. – Steve Jun 12 '17 at 09:35
  • even if b is a vector !! – Ari Jun 12 '17 at 09:41
  • `y` is 10000-by-1, `c` can be 10000-by-1, solution not valid for `x=0` – Steve Jun 12 '17 at 09:43

2 Answers2

1

Let me assume you would like to write a generic numerical solver for dy/dx = a(x)*y + b. Then you can pass the function a(x) as an argument to the right-hand side function of one of the ODE solvers. e.g.

a = @(x) 1/x;
xdomain = [1 10];
b = rand(10000,1);
y0 = ones(10000,1);
[x,y] = ode45(@(x,y,a,b)a(x)*y + b,xdomain,y0,[],a,b);
plot(x,y)

Here, I've specified the domain of x as xdomain, and the value of y at the bottom limit of x as y0.

ndalchau
  • 33
  • 3
  • I tried to simulate your code, but honestly I'm not sure I got the right result !! – Ari Jun 12 '17 at 10:21
  • Can you check the link plz – Ari Jun 12 '17 at 10:24
  • @Ar1 That looks fine, each entry of `y` has been plotted. What do you expect to see? – Steve Jun 12 '17 at 10:27
  • If for example I have this initial condition: y (x = 50) = 0 How can I add it to my code? In fact, I tried to add it but unfortunately an error msg is displayed indicating that the matrices don't have the same dimension ! Do you have any suggestion? I so am grateful because both of you helped me :D – Ari Jun 12 '17 at 13:11
  • @Ar1 you would have `xdomain = [50 100]` (or whatever your `x` domain is), `y0 = zeros(10000,1)`. – Steve Jun 13 '17 at 10:28
  • hi all here is my code load ('Bx'); load ('By'); load ('Bz'); b= sqrt(Bx.^2+By.^2+Bz.^2); fe=83333; Te=1/fe; N=length(Bz); t=0:Te:(N-1)*Te; t=t'; dB = diff(b)./diff(t); a = @(x) 1/x; xdomain = linspace(1,100,10); y0=1; for i=1:numel(dB) f=@(x,y) -(a(x)*y+dB(i)); [x,y]=ode45(f,xdomain,y0); yvec(i,:)=y(:); end figure (1) subplot(2,2,1), plot(b,'c') title('b= sqrt(Bx.^2+By.^2+Bz.^2)') subplot(2,2,2), plot(dB,'m') title('dB = -diff(b)./diff(t);') subplot(2,2,3),plot (x,y,'k') title ('solution ode45') – Ari Jun 14 '17 at 08:15
  • but i am not sure if i got the good result, here is the result of simulation – Ari Jun 14 '17 at 08:18
0

From my comments, you can solve this without MATLAB. Assuming non-zero x, you can use an integrating factor to get a 10000-by-1 solution y(x)

y_i(x) = b_i*x*ln(x) + c_i*x

with 10000-by-1 vector of constants c, where y_i(x), b_i and c_i are the i-th entries of y(x), b and c respectively. The constant vector c can be determined at some point x0 as

c_i = y_i(x0)/x_0 - b_i*ln(x0)
Steve
  • 1,579
  • 10
  • 23